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ABSTRACT 


Atmospheric  turbulence  severely  degrades  images  of  astronomical  objects. 
Providing  images  that  accurately  reflect  the  true  nature  of  these  objects  is  essential 
to  their  understanding.  Several  object  recovery  techniques  exist  within  the  field  of 
speckle  imaging  that  produce  accurate  representations  of  astronomical  objects.  This 
thesis  provides  an  in-depth  comparison  of  two  such  techniques,  Knox-Thompson  and 
triple-correlation. 

Through  computer  simulation,  this  thesis  accurately  compares  the  abilities  of 

both  recovery  techniques  to  enhance  turbulence  degraded  objects  by  exploiting  the 

di&action-limited  information  contained  in  short  exposure,  or  "speckle",  images.  The 

* 

simulation  produced  these  images  by  creating  an  object  and  several  phase  screens 
which  simulated  the  effects  of  turbulence.  Together,  the  object  and  the  appropriate 
quantity  of  phase  screens  yielded  the  required  short  exposure  images.  Application 
of  the  Knox-Thompson  and  triple-correlation  techniques  to  identical  sets  of  these 
degraded  images  produced  the  resulting  reconstructed  objects,  their  signal-to-noise 
ratios  and  their  azimuthal  RMS  phase  errors.  Comparison  of  these  three  factors  over 
several  imaging  criteria  concluded  that  the  superior  phase  recovery  technique  was 
triple-correlation. 
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I.  INTRODUCTION 


Scientific  research  of  astronomical  imagery  has  been  continual  since  the 
invention  of  the  telescope  at  the  beginning  of  the  seventeenth  century.  Initially,  the 
limiting  factor  on  astronomical  image  resolution  was  the  quality  and  size  of  the 
telescope  optics.  As  technology  progressed,  larger  telescopes  improved  to  the  point 
where  atmospheric  turbulence  (hereafter,  simpty  turbulence)  became  the  limiting 
factor  on  image  resolution. 

Turbulence  severely  limits  the  resolution  of  long  exposure  images  produced  by 
ground-based  telescope  imaging  systems.  Removal  of  turbulence  corruption  to 
improve  image  resolution  is  an  area  of  extensive  research.  Optimal  telescope  site 
location  minimized  the  effect  of  turbulence,  but  did  not  produce  the  near  difihaction- 
limited  images  desired.  Removal  of  turbulence  distortion  occurred  through  the 
development  of  statistical  methods,  called  speckle  imaging  techniques,  that  produced 
near  difffaction-limited  resolution.  The  basis  of  these  techniques  was  the  assumption 
that  turbulence  remains  essentially  stationary  during  a  short  exposure  image  of  the 
desired  object.  Though  distorted,  these  short  exposure  images  retain  diffi-action- 
limited  information. 

Determination  of  the  object’s  power  and  phase  spectra  are  separate  operations. 
Although  the  power  spectrum  can  be  directly  averaged  while  retaining  diffraction- 
limited  information,  the  phase  spectrum  cannot.  Instead,  either  the  cross-spectrum 
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(the  Knox-Thonpson  method)  or  the  bispectnim  (the  triple-correlation  method)  are 
averaged  because  difEraction-limited  information  is  retained.  Then  the  phase 
spectrum  is  recovered  from  either  the  cross-spectrum  or  the  bispentrum.  The 
combination  of  the  power  and  phase  spectra  generates  the  object’s  Fourier  spectrum 
which,  when  inverse  Fourier  transformed,  provides  the  recovered  image. 

Several  factors  affect  the  quality  of  image  reconstruction.  These  factors 
influence  both  phase  recovery  techniques  and  include  the  amount  of  turbulence,  the 
size  of  the  telescope,  and  the  light  level  of  the  object.  The  randonmess  of  turbulence 
and  the  effect  of  random  photon  noise,  make  eimct  image  reconstruction  impossible, 
however,  increasing  the  number  of  short-exptwure  images  in  the  averaging  process 
recovers  a  better  image.  This  thesis  compares  the  Knox-Thompson  and  triple¬ 
correlation  phase  recovery  techniques  under  several  different  imaging  conditions  to 
determine  their  ability  to  improve  the  signal-to-noise  ratio  (SNR)  of  a  degraded 
object. 


2 


n.  BACKGROUND 


A.  RESOLUTION 

In  the  image  reconstruction  process,  resolution  determines  image  quality. 
Resolution  is  defined  as: 

the  process  or  capability  of  making  distinguishable  the  individual  parts  of  an 
object,  closely  adjacent  optical  images,  or  sources  of  light  [Ref.  1]. 

From  a  Fourier  optics  perspective,  resolution  is  proportional  to  the  high  spatial 

fi'equency  content  of  the  imaging  system.  The  standard  for  resolution  is  based  on  the 

image  of  two  point  objects  (binary  star)  viewed  through  a  telescope.  The  image 

consists  of  two  overlapping  Airy  patterns  with  intensity 

(2,) 


where  J^(r|r)  is  the  first  order  Bessel  function  and  t|r  is  the  normalized  spatial 
frequency.  The  separation  of  the  two  first  order  fiinges  of  the  two  Airy  patterns 
constitutes  the  so-called  Rayleigh  limit  of  image  resolution 


AS  =  ^  , 


(2.2) 


(in  radians)  for  a  circular  aperture  where  X  is  the  wavelength  of  light  and  D  is  the 
telescope  diameter  [Ref.  2].  Telescope  imperfections  and  turbulence  prevent 
attainment  of  the  theoretical  limit  for  large  apertures.  Speckle  imaging  techniques 
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remove  distortions  due  to  these  factors  and  produce  images  with  resolution  near  the 
theoretical  limit. 

B.  TURBULENCE  EFFECTS 

Turbulence  produces  temporal  and  spatial  variations  in  atmospheric  density, 
temperature,  and  ind«c  of  refraction.  The  turbulence  in  the  atmosphere  perturbs  an 
image,  such  as  an  Airy  pattern  from  a  point  object  (star),  producing  an  extended 
image  referred  to  as  a  seeing  disk.  The  perturbation  randomizes  the  electromagnetic 
phase  front  of  the  object.  This  randomization  produces  angular  spreading,  image 
wandering  about  its  centroid,  and  scintillation  or  "twinkling".  Turbulence  effectively 
reduces  the  telescope’s  resolving  power  by  randomly  attenuating  the  high  spatial 
frequencies. 

C  HISTORICAL  PERSPECTIVE 

Until  roughly  1970,  attempts  at  solving  the  turbulence  distortion  problem  were 
limited  to  finding  the  ideal  telescope  site.  Generally,  the  sites  were  high  in  elevation 
and  at  locations  regarded  as  having  long  periods  of  atmospheric  stability.  Even  with 
great  care  in  site  selection,  the  typical  angular  resolution  obtained  was  approximately 
one  arcsecond,  the  maximum  resolution  attainable  with  a  12  centimeter  telescope. 
Although  phase  distortions  from  turbulence  constrained  resolution,  construction  of 
large  diameter  telescopes  provided  enhanced  light  gathering  capability. 

In  1970,  a  technique  developed  by  Labeyrie  enabled  recovery  of  near 
diffraction-limited  image  Fourier  moduli  [Ref.  3].  The  concept  that  a  long 
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exposure  image  was  blurred  by  turbulence  fluctuations  due  to  phase  spectrum 
blurring  and  not  power  spectrum  blurring,  provided  the  basis  of  Labeyrie’s  technique. 
Labeyrie  determined  that  a  short  exposure  image,  about  10  milliseconds,  would  freeze 
the  disturbance  yet  still  contain  near  diffraction-limited  information  of  the  object. 
Taking  many  short  exposure  images,  calculating  their  power  spectra,  then  averaging, 
enabled  a  diffraction-limited  estimate  of  the  object’s  power  spectrum  to  be  made. 
Labeyrie’s  technique  allowed  high  resolution  measurements  of  binary  star  separations. 
However,  lack  of  object  phase  information  prevented  faithful  image  reconstruction. 

In  1974  Knox  and  Thompson  developed  a  technique  for  retrieving  the  object’s 
phase  spectrum  [Ref.  4].  The  method  uses  the  Knox-Thompson  (KT) 
algorithm  to  provide  an  estimate  of  the  object’s  phase  spectrum  using  the  same  short 
exposure  images  required  for  the  estimate  of  the  object’s  power  spectrum.  The  KT 
method  calculates  the  average  cross-spectrum  of  the  object  in  Fourier  space  to 
determine  the  object’s  phase  spectrum.  Calculation  of  the  cross-spectrum  involves 
determining  the  average  correlation  of  spatial  frequency  pairs  displaced  from  each 
other  by  a  small  frequent^  differential.  The  average  provides  a  statistical  phase 
difference  approximation  from  which  the  object’s  phase  spectrum  is  obtained. 

In  1983  Lohmann,  Weigelt,  and  Wimitzer  developed  another  technique  for 
retrieving  the  object’s  phase  spectrum  [Ref.  5].  This  method,  referred  to  as 
triple-correlation  (TC),  also  uses  short  exposure  images  to  estimate  the  object’s  phase 
spectrum.  The  TC  method  calculates  the  average  bispectrum  of  the  object  in  Fourier 
space  to  determine  the  object’s  phase  spectrum.  Calculation  of  the  bispectrum 
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involves  determining  the  average  of  a  third  order  correlation  that  consists  of  a 
frequency  point,  a  point  shifted  by  an  offset,  and  a  difference  term.  As  with  the  KT 
technique,  the  average  provides  a  statistical  phase  difference  approximation  from 
which  the  object’s  phase  spectrum  is  obtained. 

Either  the  KT  or  the  TC  techm'que  determines  the  object’s  phase  spectrum 
which  is  necessary  to  produce  accurate  recovered  images.  Combining  the 
reconstructed  power  spectrum  and  phase  spectrum  produces  the  reconstructed 
Fourier  spectrum,  which  when  Fourier  transformed,  yields  the  recovered  image. 
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in.  THEORY 


A.  TURBULENCE  MODEL 

Both  the  KT  and  the  TC  techniques  utilize  short  exposure  images  to  recover 
the  object’s  phase  spectrum.  A  method  developed  by  Tyler  and  Fried  of  the  Optical 
Sciences  Company,  simulates  turbulence  to  produce  the  short  exposure  images 
required  to  test  these  recovery  techniques  [Ref.  6].  This  method  requires  the 
following  three  assumptions:  [Ref.  7] 

1.  Turbulence  is  represented  by  a  single  phase  screen  in  the  pupil  plane  of  the 
telescope. 

2.  Turbulence  is  isoplanatic,  that  is,  the  distortion  from  turbulence  and  the 
imaging  system  is  considered  shift  invariant  over  the  entire  image  plane. 

3.  The  images  are  quasi-monochromatic. 

With  these  assumptions,  a  single  short  exposure  image  becomes  a  convolution  in 
image  space 

i(^)  =©(.?)  ♦  s(ji)  ,  (3.1) 

where  i  (^)  is  the  short  exposure  image  intensity,  o  (x)  the  object  intensity,  and 
s  (Jl)  is  the  instantaneous  point  spread  function.  The  vector  j?  represents  the  two- 
dimensional  orthogonal  spatial  coordinates  x  and  y.  Using  the  convolution  theorem, 
equation  (3.1),  becomes  a  product  in  Fourier  space 
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J(il)  =  0(3)  •  5(3)  , 


(3.2) 


where  J(3)  is  the  Fourier  transform  of  the  short  exposure  image  intensity,  0(3) 
is  the  Fourier  transform  of  the  object  intensity,  and  5(3)  is  the  instantaneous 
incoherent  transfer  function.  The  vector  3  represents  the  two^imensional 
orthogonal  spatial  frequency  coordinates  u  and  v. 

The  point  spread  function,  s  {51) ,  and  thereby  the  incoherent  transfer  function, 
5(3) ,  represent  distortions  from  both  turbulence  and  the  imaging  system.  The 
assumption  of  stationary  turbulence  is  accmate  for  short  exposure  images. 
Consequently,  an  instantaneous  distribution  of  random  phases  (phase  screen) 
approximates  the  instantaneous  distortion  of  an  image  by  turbulence.  An  array  of 
random  numbers  filtered  by  a  power  spectral  density  function  and  corrected  for  low 
spatial  frequency  under-representation  can  simulate  this  phase  screen  [Ref.  6]. 
Therefore, 

5(3)  =  P(Xf3)€^*<‘”»  ,  (3.3) 

represents  the  instantaneous  coherent  transfer  function,  where  P(Af3)  is  the 
transfer  function  of  the  telescope,  is  the  instantaneous  turbulence  phase 

screen,  F  is  the  focal  length  of  the  telescope,  and  X  is  the  wavelength  of  light  [Ref. 
7].  Finally,  the  auto-correlation  of  the  coherent  transfer  function,  5(3) , 

5(3)  =  5(3)  ★  5(3)  .  (3.4) 

yields  the  instantaneous  incoherent  transfer  function  required  for  equation  (3.2). 
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B.  PHASE  SCREEN  PRODUCTION 

Three  common  techniques  for  producing  turbulence  phase  screens  exist.  One 
technique,  the  Fast-Fourier-transform  (FFT)  method,  creates  an  array  of  filtered 
white  noise  and  inverse  Fourier  transforms  the  array  to  real  space  providing  the 
phase  screen.  A  second  technique,  referred  to  as  the  Karhunen-Loeve  (KL) 
expansion  method,  uses  the  KL  expansion  with  a  basis  of  Zemike  polynomials  to 
represent  the  phase  screen.  The  third,  hybrid,  technique  referred  to  as  the 
Karhunen-Loeve-Fast-Fourier-Transform  (KLFFT)  method,  combines  the  best 
properties  of  both  techniques  and  manufactures  phase  screens  which  most  accurately 
represent  turbulence  distortions. 

1.  Fast-Fourier-Transform  Method 

The  FFT  method  provides  a  rapid  means  of  generating  a  phase  screen. 

#  • 

Initially,  creation  of  a  square  array  of  Gaussian-distributed  random  numbers  of  unity 
variance  provides  a  representation  of  the  phase  over  the  entire  aperture  of  the 
imaging  ^tem  being  evaluated.  The  array  amplitudes  are  filtered  in  Fourier  space 
radially  outward  from  the  origin  by  the  square  root  of  the  Kolmogorov  power  spectral 
density  function 

FJQ)  =  0.1517ro*^®  I  u  ,  (3.5) 

where  |  il  |  is  the  radial  distance  from  the  origin  in  frequency  units,  and  is  the 
coherence  diameter  [Ref  7].  The  origin  is  set  to  zero  removing  the  constant  (DC) 
term  from  the  phase  screen  before  applying  the  inverse  Fourier  transform.  Two 
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phase  screens  result  since  complex  numbers  in  the  array  consist  of  real  and  imaginaiy 
parts  which  are  entirely  distinct  and  statistically  independent. 

Though  the  FFT  method  generates  phase  screens  rapidly,  it  has 
deficiencies.  The  FFT  uses  a  finite  number  of  discrete  points,  and  consequently,  high 
and  low  spatial  firequency  cutoff  occurs.  High  spatial  frequenqr  cutoff  is  minor  since 
most  of  the  wave  firont  error  induced  by  turbulence  is  of  low  spatial  frequency.  Low 
spatial  frequent^  cutoff  is  more  serious  as  it  induces  under-representation  of  low 
spatial  frequencies  producing  wave  front  tilt,  or  centroid  position  errors.  Therefore, 
the  associated  structure  function  of  the  FFT-produced  phase  screen  does  not 
completely  represent  the  5/3  power  law  turbulence  structure  function. 

2.  Karhunen-Loeve  Expansion  Method 

The  KL  expansion  method  provides  an  accurate  method  of  generating  a 
phase  screen.  Random  phases  associated  with  turbulence  can  be  expanded  in  terms 
of  a  series  of  orthogonal  functions  ( t) , 

.  (3.6) 

The  expansion  coefficients  r^,  are  uncorrelated  Gaussian  random  variables  which 
represent  turbulence  statistics.  The  orthogonal  functions  provide  the  proper  spatial 
dependence,  thereby  allowing  the  random  phase  to  be  evaluated  anywhere  within  the 
aperture.  The  above  expansion  is  refened  to  as  the  KL  expansion.  For  a  finite  value 
of  j,  the  KL  expansion  is  the  optimum  basis  whose  eigenvalues  represent  the  energy 
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content  of  the  expansion  coefficients,  ( t) ,  and  the  total  energy  is  the  sum  of  these 
eigenvalues.  [Ref.  8]. 


Zemike  expansion  coefficients,  for  a  random  phase  screen,  are  Gaussian 
random  variables.  Unfortunately,  these  expansion  coefficients  are  correlated  and 
cannot  be  used  as  a  KL  basis  set  directly.  However,  the  Zemike  covariance  matrix 
(the  matrix  of  expansion  coefficients)  is  useful  in  determining  the  KL  expansion  for 
turbulence. 

Three  properties  justify  the  use  of  Zemike  polynomials  as  a  basis  set  for 
the  KL  expansion  to  determine  wave  front  turbulence.  Use  of  the  Zemike 
covariance  matrix  provides  the  necessary  random  variables  required  for  the  phase 
screen.  Additionally,  each  eigenvector  of  the  Z^emike  covariance  matrix  is  the 
representation  of  the  KL  function  in  terms  of  the  Zemike  polynomials.  Further,  each 
eigenvalue  of  the  2^mike  covariance  matrix  is  the  variance  associated  with  the 
corresponding  KL  expansion  coefficient.  Wave  front  error  induced  by  turbulence  is 
an  outcome  of  the  these  properties.  [Ref.  6] 

The  eigenvectors  and  the  corresponding  eigenvalues  of  the  normalized 
Zemike  covariance  matrix  are  found  which  obey  the  relation 

CSi  *  ,  (3.7) 

where  C  is  the  covariance  matrix,  is  the  i**  eigenvector  and  is  the 
corresponding  normalized  eigenvalue.  The  eigenvectors  are  normalized  so  each 
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element  of  the  eigenvector  indicates  the  amount  of  the  corresponding  Zemike 
polynomial  that  is  contained  in  the  i  ^  KL  function  expressed  as 

‘^^ipZp(p)  .  (3.8) 

p 

where  K^(p)  is  the  i“*  KL  function,  is  the  component  of  thei“ 
eigenvector,  and  Zp(p)  is  the  p'^  Zemike  polynomial.  Hence,  the  random  wave 
front  error  produced  by  turbulence,  (i') ,  is 

♦  (f)  =  ,  (3.9) 

where  Yi  is  a  set  of  Gaussian-distributed  random  numbers,  f  is  the  distance  from 
the  origin,  and  D  is  the  telescope  diameter. 

Though  the  KL  expansion  method  is  accurate,  deficiencies  exist,  and  care 
must  be  taken  in  its  use.  Calculating  wave  front  distortion  using  the  KL  method 
requires  an  enormously  large  number  of  Zemike  pofynomials  to  achieve  enough 
accuracy.  The  required  number  is  proportional  to  (P/rg)^.  Additionally,  numerical 
inaccuracies  exist  in  evaluating  Zemike  polynomials  of  high  radial  order.  Therefore, 
to  achieve  the  accuracy  desired,  avoidance  of  Zemike  polynomials  of  high  radial 
order  is  necessary. 

3.  Karhunen«Loeve>Fast>Fourier-Transform  Method 

The  KLFFT  method  combines  the  fast  computational  speed  of  the  FFT 
method  with  the  optimum  low  spatial  frequency  representation  of  the  KL  functions. 
This  technique  requires  a  phase  screen  to  be  developed  by  the  FFT  method.  The 
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first  five  KL  functions  are  produced  and  applied  to  this  phase  screen  to  compensate 
for  the  under-representation  of  low  spatial  frequencies.  With  this  compensation,  this 
combined  technique  closely  represents  atmospheric  turbulence. 

The  KLFFT  method  is  quite  powerful.  Since  only  five  KL  functions  are 
produced  and  applied,  the  total  number  of  operations  per  phase  screen  is 
approximately  double  that  of  the  FFT  method  alone.  Therefore,  this  phase  screen 
production  technique  is  relatively  fast.  [Ref.  6] 

C  IMAGE  RECOVERY 

The  image  recovery  techniques  represent  methods  for  providing  the 
reconstructed  image  from  a  turbulence-distorted  object  by  utilizing  several  short 
exposure  images  of  the  object  and  a  nearby  star.  The  Labeyrie  technique  recovers 
the  object’s  power  spectrum  [Ref.  3].  The  power  spectrum  provides  the  modulus  of 
the  Fourier  transform  of  the  object,  the  first  part  of  the  complex  quantity  required. 
Both  the  KT  and  TC  techniques  recover  the  object’s  phase  spectrum.  Inverse 
Fourier  transforming  the  product  of  the  modulus  and  the  phase  provides  a 
reconstructed  image  of  the  original  object. 

1.  Power  Spectrum  Recovery 

Equation  (3.2)  represents  the  image  intensity  for  a  single  short  exposure 
image  in  Fourier  space.  The  time  average  power  spectrum  of  several  short  exposure 
images  in  Fourier  space  is 
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<|I(0)|2>  =|0(0)|2  •  <|5(0)|2>  . 


(3.10) 


where  the  term  on  the  left  is  the  time  average  power  spectrum  of  the  image,  the  first 
term  on  the  right  is  the  object’s  power  spectrum  and  the  second  term  on  the  right  is 
the  incoherent  transfer  function.  Several  images  of  a  star  under  similar  imaging 
conditions  as  the  object  determine  this  transfer  function,  which  is  the  time  average 
of  the  instantaneous  transfer  functions.  The  object  and  the  star  need  not  be  in  the 
same  isoplanatic  patch  as  long  as  the  second  order  statistics  of  the  transfer  function 
are  the  same  for  both  sets  of  exposures  [Ref.  9].  Solving  for  the  object  power 
spectrum,  equation  (3.10)  becomes 


|0(u)|* 


^  <|J(g)|^> 
<|5(0)12> 


which  recovers  the  object’s  Fourier  modulus. 


(3.11) 


2.  Phase  Spectrum  Recoveiy 

Both  the  KT  and  TC  phase  recoveiy  techniques  recover  the  object’s  phase 
spectrum  by  using  the  cross-spectrum  and  bispectrum  averaging  processes 
respectively.  To  recover  the  object’s  phase  spectrum,  the  KT  technique  calculates  the 
average  cross-spectrum  while  the  TC  technique  calculates  the  average  bispectrum. 


a.  Technique  Analysis 

The  KT  technique  is  the  simpler  of  the  two  phase  recovery  methods. 
Determining  the  cross-spectrum  is  at  the  heart  this  algorithm.  The  cross-spectrum 
is  defined  as  the  product  of  two  quantities  in  Fourier  space:  an  array  point  and  the 
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complex  conjugate  of  an  array  point  which  is  shifted,  in  frequency,  from  the  original 
array  point  by  a  small  offset,  Au.  The  cross-spectrum,  €5(3) ,  is 

€5(3)  =  Jf(u)  •  Jf(a  +  AiJ)  ,  (312) 

where  Xi3)  is  the  array  point  and  X*  (3*  A3)  is  the  complex  conjugate  of  the  array 
point  shifted  by  the  ofrset  vector  A  3. 

The  TC  technique  is  a  more  complicated  form  of  phase  recovery.  The 
bispectrum  is  defined  as  the  product  of  -hree  quantities  in  Fourier  space:  an  array 
point,  the  complex  conjugate  of  an  array  point  which  is  shifted,  in  frequency,  from 
the  original  array  point  by  an  offset,  and  an  array  point  which  is  a  function  of  the 
offset  only.  The  bispectrum,  BS  (3) ,  is 

BSi3)  »  X(3)  •  X*(3*A3)  •  X(A3)  ,  (3.13) 


where  X(A3)  is  the  array  point  which  is  a  frinction  of  offset  only  and  the  other 
terms  are  the  same  as  defined  for  equation  (3.12). 

From  equation  (3.2),  the  average  cross-spectrum  is 


<I(3)  •  r{3*A3)>  = 

[0(3)  •  0*(a  +  A3)]  •  <5(3)  •  5*(3+A3)>  . 


(3.14) 


where  the  first  term  on  the  right  is  the  object’s  cross-spectrum  and  the  second  term 
on  the  right  is  the  average  incoherent  transfer  function  cross-spectrum.  The  average 
bispectrum  is 
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(3.15) 


<i{m  •  j*(0+Aa)  •  j(Aa)>  = 

[0(a)  •  o*(a+A3)  •  oi^a)]  • 

<5(3)  •  5*(3+Aa)  •  5(A3)>  , 

where  the  first  term  on  the  ri^t  is  the  object’s  bispectrum  and  the  second  term  on 
the  right  is  the  average  incoherent  transfer  function  bispectrum. 

The  average  incoherent  transfer  function  cross-spectrum  and 
bispectrum  contain  distortions  firom  both  the  turbulence  and  the  imaging  ^tem.  The 
average  phases  resulting  from  these  calculations  for  turbulence  are  zero. 
Imperfections  in  the  imaging  system  produce  phases  which  are  negligible  for  both  the 
cross-spectrum  [Ref.  10]  and  the  bispectrum  [Ref.  5]  calculations.  Therefore, 
the  average  cross-spectrum  and  bispectrum  of  the  incoherent  transfer  function  is 
assumed  to  be  unity.  Equations  (3.14)  and  (3.15)  reduce  to  simpler  forms 

0(3)  •  0*(a+A3)  -  <r(3)  •  r*(3+A3)>  ,  (3.16) 

and 


0(3)  •  0*(3  +  A3)  •  0(A3)  « 

(3.17) 

<J(3)  •  J*(3+A3)  •  J(A3)>  . 
b.  Phasor  Spectrum  Recovery 

Direct  recovery  of  the  object’s  phase  spectrum  is  not  possible  since 
the  cross-spectrum  and  bispectrum  phases  are  only  known  modulo  In.  The  recursion 
algorithm  fails  when  the  cross-spectrum  or  bispectrum  phase  estimates  are  not  equal 
to  their  principle  arguments  if  multiple  estimates  of  a  single  object  phase  spectrum 
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point  are  used.  To  avoid  this  problem  without  losing  information,  the  reconstructed 
object’s  phasor  spectrum  is  determined  instead.  The  recovery  process  wfll  henceforth 
be  referred  to  as  phasor  reconstruction. 

An  arbitrary  complex  number  N  in  phasor  notation  is 

^  =  I  I  ,  (3.18) 

where  |  |  is  the  modulus  of  the  complex  number  and  is  the  phasor  in  which  <)> 

represents  the  phase  of  the  complex  number.  Solving  for  the  phasor,  equation  (3.18) 
becomes 


e‘*  •  - 


N 

\N\ 


(3.19) 


where  the  subscript  ph  denotes  phasor.  Therefore,  solving  for  the  phasors  of 
equations  (3.16)  and  (3.17)  by  dividing  each  their  respective  moduli  gives 

•  0^(3+An)  =  ,  (3.20) 

and 

'  Op^(AU)  =  i^(i2,Aa)  ,  (3.21) 

where 


<j(u)  •  r{ii*Aa)> 
|<I(u)  •  r*(iI+A0)>|  ' 


(3.22) 


and 


rSca.AO)  .  <J(a)-J-(g^AO),-_j(Aa)> 

^  |<I(a)  •  J*(il+A0)  •  J(A0)  >1 


(3.23) 
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aTic  tenns  J^(Q, AO)  and  l^(U,AiD  are  four-dimensional  quantities,  since 
several  offset  values  may  be  used  to  determine  estimates  of  the  object’s  cross¬ 
spectrum  and  bispectrum.  Solving  for  the  ofiset-shifted  object  phasor  and  applying 
the  complex  conjugate  operator  to  both  sides,  equations  (3.20)  and  (3.21)  become 


0^(a*Aa) 


'  I^ia.AQ)  ] 
) 


(3.24) 


and 


o^{a*Am 


(  i^ia.Arn  V 

^o^(3)-o^(Aa)J  • 


(3.25) 


Equations  (3.24)  and  (3.25)  provide  the  phasor  spectrum  necessary  for  image 
reconstruction.  However,  use  of  these  equations  is  limited  to  infinite  photon  count, 
short  exposure  images  which  are  unrealistic  and  useful  for  computer  simulations  only. 

c.  Photon  Noise  effects 

Compensating  for  photon  noise  allows  image  reconstruction  firom  low 
photon  count  images,  though  more  of  these  images  are  required  in  the  process  to 
resolve  the  object.  During  the  recovery  process  of  low  light-level  objects,  the 
introduction  of  photon  bias  occurs.  This  bias  corresponds  to  the  cross-spectrum  or 
bispectrum  averaging  of  photon  events  with  themselves  and  contributes  no  useful 
information  to  the  average.  With  the  bias  removed,  equations  (3.11)  and  (3.12) 
become 
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cs(a)  ^  x(a)  ^  xuu+AW  -  x*(Aa)  , 


(3.26) 


and 


BS(a)  =  X(a)  'XUO  +  All)  ‘  X(A3)  +  2Np  - 
+  iJf(3+Aa)|2  _  |x(Aa)|2  , 


where  -j:*(AiI)  removes  the  cross-spectrum  photon  noise  bias, 

-\X{0  +  AQ)|^,  and  -|jr(Au)|^  remove  the  bispectrum  photon  noise  bias,  and  Np 
is  the  photon  count  [Ref.  11].  Therefore,  for  realistic  image  reconstruction, 
equations  (3.24)  and  (3.25)  become 


Op^{a*Am 


rXTBZAS 

J-ph 


(Q.Afl) 


and 


Opft(a+Aa) 


(  i^^(a,Am  y 
[op,im-o^ua)j 


(3.28) 


(3.29) 


where 


jjcTBiAS/0  _  <J(u)  •  J*(a-»-Ag)  -  J*(Ag)> 

^  '  |<J(u)  •  J*(a+Au)  -  J*{Au)>|  '  (3.30) 


and 


-TCBZAS 

•Lpb 


(a.  Am  = 


<j(a)  •  j*(a+Aa)  -  iKOip-  |j(a+Aa)i*-  \i{Am\^*  2n^> 
\<i(m  •  r(G+A0)  -  |J(G)|2-  ij(a+Au)i2-  |j(ag)|2+  2Np> 


(3.31) 
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Equations  (3.28)  and  (3.29),  allow  realistic  image  reconstruction. 


d.  Phasor  Spectrum  Wei^iting 

Phasor  spectrum  weighting  provides  a  means  for  enhancing  desired 
phasor  estimates,  increasing  reconstructed  image  quality.  Both  phasor  recovery 
techniques  benefit  from  phasor  spectrum  weighting.  Weighting  techniques  suppress 
higher  frequencies  and  enhance  the  reconstructed  image’s  SNR.  The  method 
presented  is  the  weighted  least  squares  estimation  approach.  Matson  showed  that 
this  method  obtained  the  best  results  of  four  approaches  anatyzed  [Ref.  12]. 
The  method  weights  the  object’s  phasor  spectrum  with  the  SNR  determined  from  the 
variance  of  the  cross-spectrum  or  bispectrum  as  follows: 


a^dmlSSiQ)]) 


n 

s 

l?e[55j(3)2]  -  {12e[<55(a)>]}2^ 

n-l 

n  -  1 


(3.32) 


(3.33) 


o  *  V®*  Jjn[55(u)  ] ) 


(3.34) 


sm  *  ^  , 

o 

where  55(3)  represents  either  the  cross-spectrum  or  the  bispectrum. 


(3.35) 


20 


e.  Recovery  Process 

Determination  of  the  object’s  phasor  spectrum  is  afforded  by  the 
recursive  method.  An  estimate  of  the  image  cross-spectrum  or  bispectrum  for  each 
short  exposure  image  spectrum  element  provides  the  necessary  information  to 
determine  the  object’s  phasor  spectrum.  The  estimates  for  each  array  element  are 
determined  recursively,  from  the  origin  radially  outward  to  the  diffraction  limit  of  the 
imaging  ^tem.  This  method  capitalizes  on  the  inherent  property  of  the  phasor 
spectrum  SNR  which  decreases  radially  as  a  frinction  of  increasing  spatial  frequency. 
These  estimates  are  averaged  over  all  the  short  exposure  images  and  are  weighted 
utilizing  equation  (3.3S).  This  process  results  in  the  object’s  phasor  array  determined 
by  equations  (3.28)  or  (3.29).  Once  the  reconstructed  phasor  array  is  found  in  this 
manner,  it  is  multiplied  by  the  square  root  of  the  power  spectrum  provided  by 
equation  (3.11).  Inverse  Fourier  transforming  this  spectrum  to  image  space  provides 
the  final  reconstructed  image. 

D.  PHASOR  RECOVERY  TECHNIQUE  COMPARISON 

There  are  several  distinctions  between  the  two  phasor  recoveiy  techniques, 
aside  from  the  obvious  differences  of  equations  (3.28)  and  (3.29).  Since  the  recursion 
technique  utilizes  all  the  values  from  unity  to  the  offset  value  in  the  averaging 
process,  the  greater  the  offset,  the  better  the  reconstructed  image.  For  the  KT 
technique,  the  modulus  of  the  offset  vector,  |  dd  |,  is  restricted  to  be  less  than  the 
seeing  limit,  Tq/X,  whereas  for  the  TC  technique,  the  offset  is  unlimited.  This 
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restriction  allows  the  TC  technique  to  surpass  the  KT  technique  in  image  quality  by 
using  a  greater  offset  that  includes  more  phase  information.  [Ref.  13] 

Shift  invariance  becomes  important  for  low  photon  count  images.  The  KT 
technique  is  not  shift  invariant  and  requires  shifting  of  the  degraded  short  exposure 
images  to  the  center  of  the  image  array  by, 


and 


X  = 


•  Ij(x) 
Ij(x) 


(3.36) 


yj  •  Ij  (y) 


(3.37) 


prior  to  the  recovery  process.  Image  centroiding  minimizes  linear  phase  ramping  in 
Fourier  space,  however,  centroiding  accuracy  decreases  with  lower  light  levels 
resulting  from  the  randomness  of  the  image.  Consequently,  since  the  KT  technique 
is  not  shift  invariant,  it  performs  poorly  for  low  photon  count  images  because  the 
centroiding  accuracy  is  decreased.  The  TC  technique  is  shift  invariant,  eliminating 
the  centroiding  requirement.  Since  the  TC  technique  does  not  require  centroiding, 
the  inaccuracy  of  the  centroiding  process  at  low  light  levels  does  not  enter  into  TC 
image  reconstruction.  [Ref.  11] 
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E.  SUMMARY 


The  short  exposure  image  is  the  key  to  effective  image  recovery.  The  high 
spatial  frequencies  within  the  short  exposure  image  contain  diffraction-limited 
information,  though  the  object’s  phase  spectrum  is  randomized.  Recovery  of  the 
phasor  spectrum  provides  the  information  required  for  true  image  reconstruction. 

Three  techniques  for  image  reconstruction  provide  the  recovered  image.  The 
Labeyrie  technique  recovered  the  object’s  power  spectrum.  Both  the  Knox- 
Thompson  and  the  Triple-Conelation  techniques  recovered  the  object’s  phasor 
spectrum.  Each  phasor  recovery  technique  yielded  an  independent  array  of  phasors 
that,  when  combined  with  their  corresponding  power  spectrum  and  inverse  Fourier 
transformed,  produced  reconstructed  images  for  comparison.  Simulation  of 
turbulence-degraded  short  exposure  images  by  the  Karhunen-Loeve-Fast-Fourier- 
Transform  method  allowed  computer  simulated  comparison. 
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IV.  COMPUTER  SIMULATION 


A.  COMPUTER  REQUIREMENTS 

Turbulence  simulation  and  image  recoveiy  processing,  by  their  nature,  require 
enormous  amoimts  of  calculations.  The  simulation  presented  in  this  thesis  creates 
an  object  and  a  turbulence  phase  screen.  The  phase  screen  distorts  the  object 
producing  a  short  exposure  image.  Several  of  these  short  oposure  images  are 
utilized  in  the  image  reconstruction  process  by  using  either  the  KT  or  TC  recovery 
techniques.  The  phase  screens  and  short  exposure  images  are  presented  in  the  form 
of  two-dimensional  square  arrays.  The  arrays  produced  are  of  dimension  64  x  64 
consisting  of  4096  elements.  Several  operations  are  performed  on  each  element  in 
these  arrays  throughout  the  entire  simulation  process.  As  a  result,  the  requirement 
arose  for  a  fast  computer  to  process  the  arrays  and  for  a  large  random  access 
memoiy  (RAM)  to  store  them  during  the  process. 

A  personal  computer  was  used  for  the  simulation  process  and  the  data 
reduction.  The  computer,  a  Compaq  Deskpro  80386/20  with  16  megabytes  of  RAM 
and  a  Weitek  1167  coprocessor,  provided  ample  speed  and  convenience  as  long  as 
array  sizes  did  not  exceed  64  x  64.  Standard  Fortran  77  was  the  language  used 
throughout  the  simulation.  A  Microway  32  bit  NDP  Fortran-386  compiler  provided 
the  speed,  precision,  and  array  processing  capabilities  required  for  the  simulation. 
The  construction  of  each  short  exposure  image  and  its  subsequent  cross-spectrum  or 
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bispectrum  estimation,  required  approximately  28  to  34  seconds  to  process  with  this 
computer-compfler  combination,  depending  on  the  imaging  parameters  (coherence 
length,  photon  count,  short  exposure  image  quantity).  Surfer  and  Grapher  software 
from  Golden  Software  provided  plots  of  image  data,  phasor  spectrum  SNR 
comparison  data,  and  azimuthal  RMS  phase  error  (henceforth,  simply  phase  error) 
comparison. 

B.  SIMUIATION  PROCESS 

The  simulation  process  compared  the  KT  and  TC  phasor  recovery  techniques. 
The  necessity  for  two  programs,  one  using  the  KT  technique  and  the  other  the  TC 
technique,  arose  from  RAM  limitations.  To  ensure  accuracy,  both  programs  were 
identical  accept  for  the  individual  phasor  recovery  subroutines.  Additionally,  the 
imaging  parameters  and  the  random  number  seeds  were  identical  for  each 
comparison  run  to  ensure  production  of  the  same  phase  screens  and,  hence,  the  same 
short  exposure  images.  With  identical  short  exposure  images  and  object  power 
spectra,  the  only  distinguishable  difference  between  reconstructed  images  resulted 
from  the  utilization  of  different  phasor  recovery  techniques. 

The  simulation  utilized  the  speckle  imaging  procedure.  This  procedure  involved 
the  use  of  several  short  exposures,  from  25  to  1600,  to  remove  the  effects  of 
turbulence  by  means  of  an  averaging  process.  This  process  determined  the  object’s 
power  and  phasor  spectra  by  calculating  the  autocorrelation  and  the  cross-spectrum 
or  bispectrum  respectively,  for  each  short  exposure  image.  At  the  end,  these  values 
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were  averaged  and  combined,  providing  the  object’s  Fourier  spectrum.  A  separate 
program  filtered  and  transformed  this  spectrum  to  image  space,  yielding  the 
recovered  image.  The  simulation  process  is  shown  in  Figure  4.1. 

1.  Object  Production 

Construction  of  an  object  that  provided  adequate  detail  to  test  the 
resolution  of  the  two  phasor  recovery  techniques  was  essential.  Three  objects  were 
designed  to  compare  the  two  techniques  over  several  different  image  parameters. 
These  objects  were  created,  scaled,  transformed  to  Fourier  space,  and  normalized. 

a.  Object  Creation 

The  first  step  in  the  process  created  the  object.  The  option  to 
construct  one  of  three  objects  was  provided.  The  first  object  resembled  a  finite¬ 
dimensional  astronomical  body  centered  in  the  array.  The  Ixxfy  was  a  convolution 
of  a  Gaussian  function  and  a  circular  pupil  function,  giving  it  the  appearance  of  a 
smooth  planet  Since  phasor  spectrum  SNR  declines  radially  with  spatial  frequency, 
roimd  objects  provide  a  less  than  ideal  choice  for  image  recovery  comparison. 
However,  increasing  the  detail  on  the  body  provided  a  means  to  test  the  resolution 
capabilities  of  the  two  phasor  recovery  techniques.  Seven  Gaussian  functions  of 
various  size  and  depths  at  random  locations  on  the  body  provided  craters.  These 
craters  gave  the  body  the  appearance  of  an  asteroid.  The  randomness  of  the  craters 
on  the  asteroid  provided  the  additional  detail  to  test  resolution. 
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)EGIN  SIMUIATIOt 


Figure  4.1  Simulation  Process 
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The  second  object  resembled  a  binary  star.  This  object  consisted  of  two 
delta  functions  placed  symmetricaUy  about  the  center.  Initially, this  object  provided 
the  ability  to  troubleshoot  the  program  since  the  phase  spectrum  of  an  equal  intensity 
binary  star  involves  a  square  wave  pattern  that  was  easily  recognizable.  After  the 
completion  of  troubleshooting,  the  binary  star  allowed  a  test  of  the  ability  of  the 
phasor  recovery  techniques  to  resolve  point  objects  at  very  low  photon  counts.  One 
point  had  twice  the  intensity  of  the  other  to  reduce  any  ambiguities  brought  about  by 
symmetry. 

The  third  object  resembled  a  star.  The  object  consisted  of  a  delta 
function  at  the  center  of  the  array.  The  short  exposure  image  of  a  point  source  yields 
the  instantaneous  incoherent  transfer  function  representing  distortion  from  the 
turbulence  and  the  imaging  system.  The  transfer  function  allowed  recovery  of  the 
object’s  power  spectrum. 

b.  Object  Scaling 

The  object  had  to  be  scaled  before  use.  Object  scaling  assured  that 
the  imaging  parameters  retained  complete  frequency  information  within  the  given 
array  size  and  ensured  that  the  size  of  the  object  was  within  acceptable  limits.  One 
such  parameter  was  rQ,  the  coherence  length,  which  was  a  measure  of  the  amount  of 
turbulence  present  in  the  atmosphere  [Ref.  14].  For  the  simulation,  values 
of  0.206  and  0.103  meters  were  chosen.  Another  parameter,  |u|,  was  the  offset 
value.  This  value  represented  the  number  of  pixels  (array  elements),  in  Fourier 
space,  contained  within  the  coherence  length.  The  offset  maximized  the  averaging 
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of  cross-spectrum  and  bispectrum  while  preventing  loss  of  high  frequency  information. 
O&et  values  of  two  and  four  were  chosen  for  the  simulation.  The  corresponding 
width  of  each  pixel  in  terms  of  coherence  length  was  determined  and  divided  into  the 
telescope  diameter  to  calculate  the  number  of  pixels  retained  by  this  telescope  imder 
the  conditions  of  the  above  parameters.  The  number  of  pixels  is  synonymous  with 
the  frequency  cutoff  of  the  diffraction-limited  telescope  incoherent  transfer  function 
(.D/X) .  The  image  array  size  limits  the  frequency  cutoff  value  to 


(4.1) 


where  is  the  frequency  cutoff  and  n  is  the  array  size.  If  the  frequency  cutoff  is 
too  large,  frequency  information  is  lost. 

After  constraining  the  imaging  parameters,  the  field  of  view  and  the  object 
size  were  ascertained.  The  field  of  view  (FOV)  is  equivalent  to  the  reciprocal  of  the 
fundamental  spatial  frequency 


FOV  =  |iJ|  • 


(4.2) 


where  (X/r^)  is  the  seeing  disk.  In  general,  the  allowance  of  one  seeing  disk  width 
between  the  object  and  each  side  of  the  array  provided  for  object  distortion  effects. 
The  size  of  the  asteroid  was  then  maximized  to  provide  the  most  visual  detail  while 
satisfying  the  above  criteria.  Verification  that  the  binary  star  and  the  star  obeyed  the 
above  criteria  was  sufficient  for  those  objects. 
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c.  Fourier  Spectrum  Determination 

Production  of  phase  screens  occurred  in  Fourier  space,  while  object 
production  transpired  in  image  space.  Proper  image  production  required 
transformation  to  frequency  space.  A  two-dimensional  Fast-Fourier-transform  (2D- 
FFT)  algorithm  transformed  the  object  to  Fourier  space.  The  FFT  provides  a  fast 
and  accurate  method  of  transforming  a  discrete  frinction  to  Fourier  space.  The  2D- 
FFT  employs  a  ID-FFT  provided  by  Gonzalez  and  Wintz  [Ref.  15].  This 
ID-FFT  determines  the  discrete  Fourier  transform  of  a  complex  one-dimensional 
array  of  numbers.  The  2D-FFT  simply  calls  the  ID-FFT  for  each  row  then  each 
column  of  the  two-dimensional  object  array.  Testing  of  the  2D-FFr  by  transforming 
a  normalized  pupil  function  then  inverse  transforming  enabled  comparison  between 
the  results  and  the  ori^al  function.  With  double  precision  complex  numbers,  the 
2D-FFT  provided  accuracy  to  ten  significant  figures.  In  addition  to  determining  the 
object’s  Fourier  spectrum,  this  2D-FFT  provided  the  means  for  transformation,  from 
image  space  to  Fourier  space  and  back,  extensively  throughout  the  simulation. 

d,  Fourier  Spectrum  Normalization 

Normalization  of  the  object’s  Fourier  spectrum  produced  the  correct 
number  of  photons  in  the  short  exposure  image.  In  reality,  each  short  exposure 
image  furnishes  a  photon  count,  however  for  the  simulation,  the  some  value  was 
chosen  for  all  short  exposures  used  for  each  image  reconstruction  run.  Dividing  the 
object’s  Fourier  spectrum  by  the  photon  count  normalized  the  spectrum.  Since  the 
phase  map  was  normalized  by  its  DC  value,  the  photon  count  was  the  same  for  the 
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object  and  for  the  short  exposure  image.  The  normalization  was  impoitant  for 
Poisson  noise  generation  introduced  later  in  the  process. 


2.  Turbulence  Phase  Screen  Production 

Construction  of  a  turbulence  phase  screen  that  correctly  resembled  true 
atmospheric  turbulence  allowed  accurate  phasor  recovery  technique  comparison. 
Testing  the  KLFFT  method  of  phase  screen  production  showed  it  represented  the 
actual  5/3  power  law  structure  function  closely  [Ref.  16].  Therefore,  this 
method  was  adapted  to  produce  the  phase  screens  in  this  simulation. 

a.  Gaussian  Random  Number  Array 

Each  phase  screen  involved  an  array  of  random  numbers.  The 
random  numbers  represented  the  random  phases  produced  by  turbulence.  A 
Gaussian  distributed  random  number  generator,  provided  by  the  subroutine  Gauss, 
produced  the  required  random  numbers  that  represented  the  randomness  of 
turbulence  statistics. 

b.  Filter  Function 

The  simulation  used  a  filter  function  that  represented  the  square  root 
of  the  Kolmogorov  power  spectral  density  function,  equation  (3.5).  This  function 
represented  turbulence  statistics  and  filtered  each  array  of  Gaussian  distributed 
random  numbers  to  provide  the  turbulence  structure  function.  The  resulting  array 
elements  when  put  in  terms  of  the  Rytov  approximation,  modelled  the  5/3 
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power  law  structure  function,  and  hence  provided  a  model  for  turbulence,  though  the 
low  spatial  frequencies  were  under-represented. 

c.  Karhimen-Loeve  Functions 

The  simulation  iised  the  first  five  KL  functions.  Each  filtered  array 
was  inverse  Fourier  transformed  and  then  the  five  KL  functions  were  applied  in 
image  space  to  compensate  for  the  low  spatial  frequency  under-representation.  This 
application  involved  more  than  simple  multiplication.  The  inner  product  of  the  KL 

* 

functions  and  the  filtered  array  provided  scalars  which  expressed  the  amount  of  each 

individual  KL  function  contained  in  the  array.  These  amounts  were  subtracted  from 

the  array.  The  technique  involved  random  numbers  with  variances  equal  to  the 

eigenvalues  associated  with  their  respective  KL  function.  Multiplying  the  random 

numbers  by  their  corresponding  KL  functions  and  adding  the  result  to  the  array  gave 

« 

the  corrected  phase  screen.  The  resulting  array  elements  when  put  in  terms 

of  the  Rytov  approximation,  accurately  depicted  the  5/3  power  law  structure  function. 

An  effect  of  the  KLFFT  method  was  the  inclusion  of  tilt  in  the  short 
exposure  image.  An  option  was  given  in  the  simulation  which  allowed  the  removal 
of  tilt  setting  the  first  two  KL  functions  to  zero.  This  option  allowed  phasor 
recovery  without  the  presence  of  tilt  and  was  utilized  as  a  criterion  for  comparison  * 

of  the  phase  recovery  techniques. 
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d.  Incoherent  Tranter  Function 

Use  of  the  Rytov  approximation  determined  the  incoherent  transfer 
function  of  the  turbulence  and  imaging  system.  The  Rytov  approximation, 

,  (4.3) 

provides  the  coherent  transfer  function  of  the  turbulence  and  imaging  system,  where  A 
represents  the  amplitude  (set  to  one)  and  represents  the  realization  of  the  phase 
screen  determined  above.  Calculating  the  autocorrelation  of  the  coherent  transfer 
function  produced  the  incoherent  transfer  function.  This  calculation  first  required 
multipfying  the  phase  screen  array  produced  by  the  Rytov  approximation  by  the  pupil 
function  of  the  telescope.  Squaring  the  modulus  of  the  Fourier  transform  of  this 
array  and  Fourier  transforming  and  normalizing  the  result  yielded  the  incoherent 
transfer  function  of  the  true  phase  screen. 

3.  Object  Degradation 

The  product  of  the  phase  screen  and  the  object’s  Fourier  spectnun  yielded 
the  short  exposure  image  spectrum.  The  inverse  Fourier  transform  of  this  spectrum 
resulted  in  the  required  short  exposure  image.  Since  the  object  was  normalized  to 
the  photon  count  of  the  short  exposure  image  and  the  phase  screen  was  normalized 
to  unity,  the  final  step  in  the  degradation  process  was  to  apply  photon  noise  to  the 
short  exposure  image.  The  photon  noise  effect  was  added  to  the  short  exposure 
image  by  entering  each  image  element  into  a  Poisson  distributed  random  function 
generator  provided  by  the  Poisson  function  in  the  simulation.  The  generator  returned 
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a  random  number  at  each  image  point  drawn  from  a  Poisson  distribution  with  mean 
equal  to  that  value.  As  the  photon  count  decreased  from  infinity,  the  randomness  of 
the  returned  values  increased  which  caused  the  grainy,  photon  noise  effect  With 
photon  noise  included,  the  short  exposure  image  was  complete. 

4.  Centroiding 

The  simulation  supported  centroiding  the  degraded  short  exposure  in 
image  space.  If  desired,  the  image  array  was  shifted  to  its  true  centroid  first  by 
column,  then  row,  based  on  the  centroids  determined  from  equations  (3.32)  and 
(3.33).  Determining  the  centroid  of  the  binary  star  object  modified  with  equal 
intensity  stars  checked  the  accuracy  of  the  centroid  subroutine.  The  object  initially 
had  an  arbitrary  translation  from  its  centroid.  Centroiding  this  translated  object  using 
the  subroutine  then  analyzing  its  phase  spectrum  ensured  the  subroutine  performed 
correctly.  Image  centroiding  worked  well  on  high  light-level  images  of  the  asteroid 
above  10^  photons.  Images  below  this  level  had  unevenly  distributed  intensities  and 
centroiding  was  ineffective. 

5.  Image  Recovery 

The  simulation  separated  image  recovery  into  two  distinct  parts.  First,  the 
Labeyrie  technique  provided  the  object’s  power  spectnun.  Use  of  this  technique  for 
both  programs  ensured  uniformity  in  comparison.  Second,  the  KT  and  TC  techniques 
reconstructed  the  object’s  phasor  spectrum,  each  in  separate  programs. 
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a.  Power  Spectrum  Recovery 

Recovery  of  the  object’s  power  spectrum  provided  the  modulus  of  the 
object’s  Fourier  spectrum.  Estimates  of  the  autocorrelation  were  calculated  for  each 
short  exposure  image  of  the  object  and  the  point  source.  Averages  of  these  results 
yielded  the  average  power  spectrum.  Normalizing  both  power  spectrum  arrays  by 
their  respective  DC  values  provided  an  equivalent  intensity  basis.  Dividing  the 
object’s  average  power  spectrum  by  the  point  source’s  average  power  spectrum 
removed  imaging  system  errors.  The  square  root  of  this  result  multiplied  by  the 
photon  count  furnished  the  object’s  modulus. 

b,  Phasor  Spectrum  Recovery 

The  object’s  phasor  spectrum,  with  its  modulus,  determined  the 
object’s  Fourier  spectrum.  Each  short  e?qx}sure  image  provides  estimates  of  the 
either  the  cross-spectrum  or  the  bispectrum,  for  various  offset  values.  Several  of 
these  estimates,  each  with  a  different  offset  value,  determined  a  specific  phasor  array 
element  by  averaging  these  estimates  over  all  short  exposure  images,  then  over  aU 
offset  values. 

Recursive  calculation,  outward  fi-om  the  origin  of  the  image  array, 
supplied  the  estimates,  and  hence  the  object’s  phasor  spectrum.  The  number  of 
estimates  which  existed  within  an  estimation  circle  of  integer  pixel  radius  equal  to  the 
whole  part  of  the  pixel  distance  of  the  desired  point  determined  the  maximum 
number  of  estimates  and  the  offset  value  for  each  phasor  estimated.  The  estimation 
circle  began  at  the  origin,  and  the  estimates  of  the  phasor  spectrum  points  began  one 
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radial  pixel  distance  beyond  that.  The  phasor  at  the  origin  had  a  value  of  one,  since 
it  had  no  imaginaiy  part,  and  it  was  normalized.  Each  phasor  spectrum  point 
included  more  estimations  as  the  recursive  process  proceeded,  out  to  the  maximum 
radial  o&et  value.  The  same  estimates  of  the  phasor  spectrum  were  found  for  every 
short  exposure  image  and  then  averaged.  Averaging  a  set  of  estimates  associated 
with  a  phasor  spectrum  point  produced  the  corresponding  phasor  value  for  that  point. 

Using  the  phasor  spectrum  of  the  uncorrupted  binary  star  object 
modified  with  equal  intensity  stars  verified  the  phasor  recovery  process.  Comparing 
the  phasor  spectrum  of  the  binary  star  before  and  after  the  recovery  process  using 
complex  double  precision  numbers  provided  a  match  for  all  points  to  ten  significant 
figures. 

The  phasor  recovery  process  required  an  extensive  amount  of 
calculation.  The  time  for  phase  reconstruction  was  approximately  30  seconds  for 
every  short  exposure  image.  The  process  was  made  less  time  consuming  by  invoking 
Hermitian  symmetry.  Hermitian  symmetry  dictates 

Oi.j  =  Oli  ,  (4.4) 

therein  allowing  half  the  number  of  calculations  to  determine  the  full  object’s  phasor 
spectrum. 

At  the  end  of  the  phasor  recovery  operation,  the  variance  of  the  cross¬ 
spectrum  or  bispectrum  estimates  provided  an  SNR  value  for  each  phasor  element 
from  equations  (3.34)  through  (3.37).  The  square  of  the  SNR  for  each  estimate  was 
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then  multiplied  1^  the  corresponding  estimate  when  the  phasor  spectrum  points  were 
calculated.  This  calculation  provided  a  weighted  least  squares  estimation  of  the 
object’s  phasor  spectrum.  When  combined  with  the  recovered  modulus,  the  object’s 
Fomier  spectrum  resulted. 

6.  Azimuthal  Signal-To-Noise  Ratio 

Averaging  the  SNR  values  of  the  cross-spectrum  and  bispectrum  provided 
an  average  SNR  value  for  each  phasor  element.  This  average  SNR  array  was  then 
averaged  azimuthally,  one  radial  pixel  value  at  a  time,  from  the  origin  out  to  the 
cutoff  frequency  to  provide  a  SNR  as  a  function  of  radial  spatial  frequency.  This 
radial  SNR  provided  one  means  to  compare  the  two  phasor  recovery  techniques  as 
well  as  to  determine  the  frequency  at  which  noise  overcame  signal  to  enable  proper 
filtering. 

7.  Fourier  Spectrum  Filtering 

The  final  result  of  the  simulation  process  was  a  weighted  least  squares 
estimate  of  the  object’s  Fourier  spectrum.  Before  inverse  Fourier  transformation,  the 
object’s  Fourier  spectrum  required  filtering.  A  simple  rectangular  low-pass  filtering 
method,  which  truncated  spatizil  frequencies  beyond  the  radial  frequency  where  the 
azimuthal  SNR  was  unity,  determined  the  object’s  Fourier  spectrum.  This  filtering 
process  was  provided  by  a  separate  program  which  included  a  method  to  determine 
the  phase  error  of  the  recovered  image’s  Fourier  spectrum. 
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8.  Azimuthal  RMS  Phase  Error 

Measuring  the  phase  error  of  the  two  phasor  recovery  techniques 
produced  another  means  for  comparison.  The  object’s  phasor  spectrum  determined 
its  phase  spectrum.  Computing  the  phase  spectrum  modulo  Iv,  then  subtracting  the 
result  from  the  object  phase  spectrum  of  the  true  object  representation  provided  the 
array  point  phase  error.  The  square  of  this  error  was  determined  then  averaged 
azimuthally  one  radial  pixel  value  at  a  time,  from  the  origin  out  to  the  cutoff 
frequent^.  The  square  root  of  this  average  provided  the  azimuthal  RMS  phase  error. 

C  SUMMARY 

The  simulation  process  obtained  the  reconstructed  image  from  several  short 
mqx)sures  images.  The  process  p^odr^d  the  desired  object  and  the  phase  screens 
to  make  the  short  exposure  images  required  for  speckle  imaging.  With  several  short 
exposure  images  of  the  object,  the  Labeyrie  technique  recovered  its  power  spectrum 
and  the  Knox-Thompson  and  Triple-Correlation  techniques  recovered  its  phasor 
spectrum.  Low-pass  filtering  removed  noise  and  the  SNR  and  phase  error 
calculations  presented  means  for  comparison  of  the  two  phasor  recovery  techniques. 
The  inverse  Fourier  transform  of  the  filtered  spectrum  yielded  the  recovered  image 
presented  in  the  form  of  a  contour  plot. 
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V.  SIMUIATION  RESULTS 


A.  RECOVERY  TECHNIQUE  COMPARISON  CRITERIA 

Seven  criteria  provided  a  basis  for  comparison  of  the  KT  and  TC  phasor 
recovery  techniques.  Each  criterion  tested  the  reconstructed  image’s  resolution 
produced  by  both  algorithms  over  a  range  of  values  for  a  specific  imaging  parameter. 
A  baseline  of  imaging  parameters  (Table  A.I)  was  established.  Typically  one 
parameter  was  varied  within  an  individual  criterion.  Each  criterion  used  a  range  of 
two  to  four  imaging  parameter  values.  The  criteria  included  reconstructed  image 
evaluation  based  on  the  quantity  of  short  exposure  images,  the  short  exposure  image 
photon  count,  and  the  amount  of  turbulence.  Additionally,  the  effects  of  short 
exposure  image  tilt  on  reconstruction  as  well  as  centroiding  in  image  space  to  remove 
it,  were  weighed.  Further,  the  effect  of  offset  value  on  cross-spectrum  and 
bispectrum  estimates  and  object  size  on  image  resolution,  were  rated. 

The  techniques  were  judged  in  terms  of  image  resolution,  phasor  spectrum 
SNR,  and  phase  error.  Each  comparison  included  both  the  KT  and  TC  image 
reconstructions.  The  evaluation  included  two-dimensional  graphs  of  the  SNR  values 
for  both  the  KT  and  TC  image  reconstruction  as  well  as  their  phase  error  values. 
The  reconstructed  images  appear  with  normalized  intensities  on  two-dimensional 
contour  plots  with  hachure  marks  indicating  the  direction  of  minima.  The  SNR  and 
phase  error  values  for  each  comparison  were  plotted  against  spatial  frequent^.  The 
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maximum  frequency  value  was  the  last  point  at  which  the  SNR  had  a  value  of  one 
or  greater.  The  reconstructed  image  plots  were  visually  compared  to  the  object’s  true 
image.  By  contrast,  their  SNR  and  phase  error  values  were  compared  to  each  other. 

B.  RECOVERY  TECHNIQUE  COMPARISONS 

Appendix  A  maintains  the  results  of  the  recoveiy  technique  comparisons  within 
the  seven  criteria.  Table  A.I  delineates  the  imaging  parameters  involved  in  the 
simulation  process  and  whether  these  values  are  variable  over  these  criteria.  Figures 
A.1  and  A.2  represent  the  true  representations  of  the  asteroid  and  binary  star 
respectively.  These  figures  show  the  results  of  the  actual  computer  generated  objects 
without  turbulence  corruption,  filtering,  or  modification  resulting  from  the  imaging 
system  including  aperture  effects,  which  all  other  figures  include.  They  were  the 
reference  figures  for  the  reconstructed  images  and  the  phase  error  calculations. 
Figures  A.3  through  A22  are  plots  of  single  short  exposure  images  of  the  asteroid, 
the  binary  star  and  the  star  that  show  the  effects  of  varying  coherence  lengths  and 
photon  counts.  They  show  the  level  of  distortion  the  objects  realize  in  the  imaging 
process.  Figures  A23  through  A25  are  plots  of  long  exposure  images  that  consisted 
of  an  average  of  100  short  exposure  images  with  tilt.  The  inclusion  of  these  images 
provides  an  appreciation  for  the  necessity  of  image  reconstruction  to  acquire  an 
image  that  more  closely  resembles  the  truth. 
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1.  Short  Exposure  Image  Quantity 

Image  resolution  increases  with  the  quantity  of  short  exposure  images  used 
in  the  reconstruction  process.  Varying  the  quantity  of  images  with  no  tilt  present 
provided  the  first  comparison  criterion  for  the  recovery  techniques.  For  the  fom 
comparisons,  the  values  of  Nf  were  400  (baseline),  25, 100,  and  1600  short  exposure 
images.  All  other  baseline  parameters  remained  unchanged.  The  relevant  plots  and 
graphs  are  Figures  A.26  through  A.41.  In  all  cases,  the  visual  image  quality  of  the 
TC  recovered  images  was  superior  to  those  recovered  by  the  KT  process.  This  image 
quality  distinction  was  especially  noticeable  at  the  lower  Nf  values  of  25  and  100.  As 
the  value  of  Nf  increased,  the  distinction  decreased  to  the  point  where  it  was  only 
slightly  noticeable  at  the  Nf  of  1600.  Anatysis  of  the  phasor  spectrum  SNR  graphs 
showed,  in  general,  that  the  TC  SNR  curves  were  ofiEset  toward  approximately  ten  to 
20  percent  greater  SNR  values  than  those  of  the  KT  SNR  curves.  Further,  the  phase 
error  graphs  showed  the  error  ciirves  resulting  from  the  TC  method  were  offset 
toward  approximately  five  to  15  percent  lesser  error  than  those  of  the  KT  method. 
These  two  series  of  graphs  confirmed  that  in  all  cases,  especially  at  low  Nf,  the  TC 
technique  outperformed  the  KT  technique. 

2.  Photon  Count 

Image  resolution  increases  with  the  amount  of  short  exposure  image 
photons  present.  Varying  the  quantity  of  photons  in  the  short  exposure  images  with 
no  tilt  present  provided  the  second  comparison  criterion  for  the  recovery  techniques. 
For  the  four  comparisons,  the  photon  count  values  were  10^  (baseline),  10^,  10^,  and 
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10^  photons  where  all  other  baseline  parameters  remained  unchanged.  Figures  A.26 
through  A.29  and  A.42  through  A.S3  are  applicable  for  this  criterion.  The  visual 
quality  of  the  TC  and  KT  reconstructed  images  for  photon  counts  of  10^  and  10^  was 
almost  identical.  For  photon  counts  of  10^  and  10^,  the  KT  technique  produced 
slightly  better  reconstructed  images.  This  photon  count  distinction  of  recovery 
technique  performance  was  confirmed  by  both  the  SNR  and  phase  error  graphs.  For 
low  photon  count  short  «posure  images,  the  KT  technique  SNR  curves  were  ofEset 
toward  approximately  ten  to  20  percent  greater  SNR  values  than  those  of  the  TC 
technique,  and  the  phase  error  curves  were  general^  offset  toward  five  to  25  percent 
lesser  error  values.  Therefore,  for  low  photon  count  image  recovery  without  tilt,  the 
KT  phasor  recovery  technique  provides  slightfy  better  resolution. 

3.  Turbulence  Magnitude  and  Offset  Value 

As  the  amount  of  turbulence  increases,  the  resolution  of  the  reconstructed 
image  decreases.  Short  exposure  images  with  no  tilt  of  coherence  lengths  0.103  and 
0.206  meters  (baseline),  provided  the  third  and  fourth  comparison  criteria  for  the 
recovery  techniques.  For  the  turbulence  magnitude  criterion,  the  0.103  meter 
coherence  length  images  required  an  offset  value  of  two,  to  ensure  inclusion  of  all 
spatial  frequencies.  For  consistency,  the  0.206  meter  coherence  length  images  used 
the  same  value.  Comparison  between  offset  values  of  two  and  four  satisfied  the 
offset  criterion  for  TC  and  KT  reconstructed  images.  Off^t  values  effect  the 
quantity  of  cross-spectrum  or  bispectrum  averaging  and  the  reconstructed  image 
quality  increases  with  increasing  offset  value.  All  other  baseline  parameters  remained 


42 


unchanged.  Figures  A.54  through  A.61  showed  the  comparisons  of  the  turbulence 
criterion  and  Figures  A.26  through  A-29  and  A.58  through  A61  showed  the 
comparisons  of  the  offset  criterion.  For  the  coherence  length  case,  the  KT  technique 
provided  a  slightly  better  image  with  greater  turbulence.  The  TC  SNR  curve  was 
offset  toward  approximately  five  to  ten  percent  greater  SNR  values  relative  to  the  KT 
curve.  However,  the  KT  phase  error  curve  was  offset  toward  approximately  five 
percent  lesser  error  values  providing  the  more  resolved  image.  For  the  offset  value 
case,  the  TC  technique  produced  a  slightly  better  image,  though  the  SNR  and  phase 
error  curves  were  almost  coincident.  This  apparent  contradiction  arose  from  TC 
techniques  having  more  fi-equent^  values  above  the  unity  SNR  value,  thereby 
providing  higher  frequencies  for  the  filtering  process. 

4.  Tilt 

The  resolution  of  the  reconstructed  image  declines  with  the  inclusion  of 
tilt  in  the  short  exposure  images.  The  previous  criterion  comparisons  were  conducted 
without  the  presence  of  tilt.  The  addition  of  tilt  provides  a  more  realistic  comparison 
of  the  phasor  recovery  techniques  as  tilt  is  always  present  in  true  short  exposure 
images.  Varying  the  photon  count  of  the  short  exposure  images  with  tilt  present 
provided  the  fifth  criterion  for  comparison  of  the  recovery  techniques.  For  the  three 
comparisons,  the  photon  count  values  were  10^,  10^,  and  10^  photons  and  all  other 
baseline  parameters  remained  unchanged.  Figures  A.62  through  A.77  were 
applicable  for  this  criterion.  In  all  photon  count  cases,  the  TC  recovered  image  were 
superior.  The  KT  SNR  curves  were  generalty  offset  toward  30  to  50  percent  lesser 
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SNR  values  and  dropped  below  the  unity  value  at  lower  spatial  frequencies  than  the 
TC  curves.  Hence,  fewer  high  frequency  values  were  included  in  the  filtering  process, 
producing  a  poorer  resolution.  The  phase  error  curves  for  the  TC  recovered  images 
were  ofEset  toward  approximately  20  to  SO  percent  lesser  error  values  than  for  those 
recovered  using  the  KT  process  except  at  the  10^  photon  count.  At  this  value, 
however,  the  KT  SNR  curve  was  offret  toward  much  lower  SNR  values. 
Consequently,  the  TC  techm'que  was  found  to  be  superior  when  tilt  was  included  in 
the  short  exposure  images. 

5.  Centroiding 

Centroiding  the  short  exposure  images  prior  to  image  reconstruction 
enhances  both  the  TC  and  KT  recovery  techniques  for  high  photon  counts.  Short 
exposure  image  centroiding  provided  the  sixth  criterion  for  comparison  of  the 
recovery  techniques  by  again  varying  the  photon  count  of  the  short  exposure  images. 
With  tilt  present,  the  images  were  centroided  prior  to  reconstruction.  For  the  three 
comparisons,  the  photon  count  values  were  10^,  10^,  and  10^  photons  and  all  other 
baseline  parameters  remained  unchanged.  Additionally,  a  comparison  with  10^ 
photons  and  1600  short  exposures  with  and  without  tUt  tested  the  effects  of 
centroiding  at  high  Nf  values.  Figures  A78  through  A93  were  relevant  for  this 
criterion.  Centroiding  offered  a  two  to  five  percent  improvement  to  both  recovery 
methods  with  higher  photon  count.  At  low  photon  counts  such  as  10^  photons, 
centroiding  actually  offset  the  phase  error  curves  toward  higher  error  values.  For  the 
TC  method,  centroiding  should  have  no  effect  or  the  effect  should  disappear  with 
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large  enough  Nf  values  since  the  method  is  shift  invariant.  However,  at  10*  photons 
and  with  1600  short  exposure  images,  a  minor  improvement  in  phase  error  occurred 
with  centroiding.  Centroiding  did  not  improve  the  reconstructed  image  to  thr*  point 
of  those  recovered  images  having  no  tilt,  nor  did  centroiding  bring  the  KT  method 
results  in  line  with  that  of  the  TC  method.  However,  a  minor  improvement  in  phase 
error  occurred  for  both  methods. 

6.  Point  Objects 

The  resolution  of  the  reconstructed  image  depends  upon  its  size  and 
detail.  Resolution  of  an  object  such  as  a  binary  star  occurs  more  easily  because  of 
the  requirement  for  less  short  exposure  image  photons.  A  binary  star  with  one  star 
having  twice  the  intensity  as  its  counterpart,  provided  the  seventh  comparison 
criterion  for  the  phasor  recovery  techniques.  The  binary  star  was  corrupted  by 
turbulence  of  coherence  length  0.103  and  0.206  meters  and  with  short  exposure  image 
photon  counts  of  10^  and  10*  photons.  Figures  A.94  through  A.  109  are  germane. 
The  two  techniques  produced  identical  results  at  the  higher  photon  count  and  lower 
turbulence  values.  With  lower  photon  count  and  greater  turbulence,  the  TC 
technique  provided  greater  resolution  of  the  stars.  Both  the  SNR  and  phase  error 
curves  supported  this  fact. 

7.  Recovery  Technique  Comparison  Findings 

For  the  eventual  use  of  real  objects  in  future  applications  of  the  recovery 
techniques  where  tilt  is  inherent  in  short  exposure  images,  triple-correlation  was  the 
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superior  reconstruction  technique.  The  triple-correlation  technique  exceeded  the 
Knox-Thompson  technique  by  far  in  the  comparisons  where  tilt  was  a  concern.  The 
shift  invariance  of  the  triple-correlation  technique  provided  the  ability  to  resolve  low 
light-level  objects  where  the  Knox-Thompson  technique  failed.  Though  the  Knox- 
Thompson  technique  required  eight  percent  less  time  for  image  reconstruction,  the 
resolution  improvement  obtained  by  the  TC  approach  outweighed  the  computational 
time  efficiency  of  the  competing  technique. 
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VI.  CONCLUSION 


A.  OVERVIEW 

This  thesis  compared  the  Knox-Thompson  and  triple-correlation  phasor 
recovery  techniques.  Since  speckle  imaging  involves  extensive  calculations, 
comparison  of  the  two  techniques  required  a  powerful  computer.  The  simulation 
produced  an  object  and  a  phase  screen  for  each  short  exposure  image.  The 
diffraction-limited  information  in  these  images  allowed  reconstruction  of  the  object’s 
power  and  phasor  spectra.  Combining  these  spectra  produced  the  reconstructed 
image.  Image  reconstruction  comparison  under  seven  imaging  criteria  permitted  the 
ability  to  determine  the  superior  technique.  The  triple-correlation  technique  provided 
the  best  overaU  image  resolution.  This  judgement  stems  from  its  superiority  with 
regard  to  realistic  short  exposure  images  which  included  tilt. 

B.  OPTIMUM  IMAGE  RECOVERY  APPROACH 

This  thesis  found  that,  of  the  two  phasor  recovery  techniques  compared,  the 
triple-correlation  technique  was  the  optimum  approach  for  real  short  exposure  image 
recovery.  From  the  shift  invariance  of  the  triple-correlation  technique,  attainment 
of  20  to  50  percent  less  azimuthal  phase  error  values  occurred  when  compared  with 
the  Knox-Thompson  technique.  As  a  result,  use  of  the  triple-correlation  phasor 
recovery  technique  is  essential.  Specifically,  removal  of  the  noise  bias  compensates 
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for  the  random  photon  noise  that  is  intrinsic  to  real  images.  The  use  of  phasor  vice 
phase  recovery  is  key  to  avoid  the  phase  2ir  wrap-around  problem.  The  weighting 
of  the  phasors  determined  from  their  cross-spectrum  or  bispectrum  estimates  by  the 
least  squares  estimation  approach,  is  critical.  The  recursion  method  used  to  provide 
the  phasors  is  not  imperative,  and  other  techniques  may  be  used,  such  as  the  method 
of  least  squares  or  the  method  of  steepest  decent,  not  discussed  in  this  thesis.  With 
regard  to  imaging  parameters,  the  triple-correlation  technique  allowed  larger 
maximum  o&et  values  than  the  Knox-Thompson  technique;  D/X  instead  of  Xg/X. 
This  provides  maximum  estimate  averaging.  Based  on  the  results  with  centroiding, 
it  is  helpful  for  relatively  high  light  level  short  exposure  images.  Finally,  peak 
recovered  image  resolution  requires  the  maximum  amount  of  short  raposure  images 
practicable. 

C  FURTHER  STUDY 

This  thesis  provided  simulated  results  for  comparison  of  the  two  recoveiy 
techniques.  Application  of  the  triple-correlation  technique  to  actual,  turbulence- 
degraded  images  provides  an  avenue  of  research.  Development  and  testing  of  other 
methods  of  extracting  the  phasor  spectrum  beyond  the  recursion  method  demands 
analysis.  Reconstructed  Fourier  spectrum  filtering  processes  beyond  the  simple 
rectangular  low  pass  filtering  approach  used  herein  require  exploration. 


48 


LIST  OF  REFERENCES 


1.  Websters  Ninth  New  Collepate  Dictionary,  p.  1(X)4,  Merriam-Webster  Inc., 

1986. 

2.  Hecht,  E.,  Optics.  2nd  ed.,  pp.  422-423,  Addison-Wesl^  Publishing  Co.,  Inc., 

1987. 

3.  Lab^rie,  A.,  "Attainment  of  Difflaction  Limited  Resolution  in  Large 
Telescopes  by  Fourier  Analyzing  Speckle  Patterns  in  Star  Images,”  Astron. 
AstroDhvs..  V.  6,  pp.  85-87,  1970. 

4.  Knox,  K.  T.,  and  Thompson,  B.  J.,  "Recoveiy  of  Images  from  Atmospherically 
Degraded  Short-Exposure  Photographs."  Astrophv.  J..  v.  193,  pp.  L45-L48, 1974. 

5.  Lohmann,  A  W.,  Weigelt,  G.  P.,  and  Wimitzer,  B.,  "Speckle  Masking  in 
Astronomy:  Triple-Correlation  Theory  and  Applications"  Appl.  Optics,  v.  22, 
pp.  4028-4037,  1983. 

6.  The  Optical  Sciences  Company  Report  TR-514,  "Wave  Front  Sensing:  A  New 
Approach  to  Wave  Front  Reconstruction,"  by  G.  A  Tyler  and  D.  L.  Fried,  pp. 
87-90,  May  1983. 

7.  Weapons  Lab  Technical  Report  89-46,  "An  Atmospheric  Computer  Simulation 
for  Speckle  Imaging  Agorithm  Analysis,"  by  C.  L.  Matson  and  I.  E.  Drunzer, 
pp.  14-29,  1989. 

8.  Whalen,  A  D.,  Detection  of  Signals  in  Noise,  pp.  285-290,  Academic  Press,  Inc., 
1971. 

9.  Knox,  K.  T.,  "Image  Retrieval  from  Astronomical  Speckle  Patterns,"  J.  Opt.  Soc. 
Am..  V.  66,  no.  11,  pp.  1236-1239,  1976. 

10.  Knox,  K.  T.,  Diffraction-Limited  Imaging  with  Astronomical  Telescopes.  Ph.D 
Dissertation,  Institute  of  Optics,  University  of  Rochester,  New  York,  1975. 

11.  Ayers,  M.  J.,  Northcott,  M.  J.,  and  Dainty,  J.  C.,  "Knox-Thompson  and  Triple- 
Correlation  Imaging  through  Atmospheric  Turbulence,"  J.  Opt.  Soc.  Am.,  v.  5, 
no.  7,  pp.  963-985,  1988. 


49 


12.  Weapons  Lab  Technical  Report  1351-54,  "Weighted  Least  Squares  Phase 
Reconstruction  from  the  Bispectrum,"  by  C.  L.  Matson,  pp.  1-27, 1990. 

13.  Beletic,  J.  W.,  "Comparison  of  Knox-Thompson  and  Bispectrum  Algorithms  for 
Reconstructing  Phase  of  Complex  Extended  Objects,"  paper  presented  at  the 
ESO  conference  on  High  Resolution  Imaging  by  teterferometiy,  Garching, 
FRG,  March,  1988,  pp.  357-362. 

14.  Walters,  D.  L.,  Favier,  D.  L.,  and  Hines,  J.  R.,  "Vertical  Path  Atmospheric  MTF 
Measurements,"  J.  Opt.  Soc.  Am.,  v.  69,  pp.  828-837,  1979. 

15.  Gonzalez,  R.  C,  and  Wintz,  P.,  Digital  Image  Processing,  p.  108,  Addison- 
Wesley  Publishing  Co.,  Inc.,  1987. 

16.  The  Optical  Sciences  Compai^  Report  TR-663,  Phase  Screen  Production, 

G.  Cochran,  pp  6-8,  September  19^. 


50 


APPENDIX  A.  PLOTS  AND  GRAPHS 


The  following  table,  plots,  and  graphs  were  referenced  in  the  text. 
Figure  AJ  Imaging  Parameters. 


Imaging  Parameter  Parameter  Parameter  Parameter 

Abbrev.  Baseline  Value 

Short  Exposure  Nf  400  Variable 

Image  Quantity 

Coherence  Length  rg  0.206  (m)  Variable 

Short  Exposure  Np  10*  Variable 

Image  Photon 

Quantity 

Offset  Value  OS  4  (pixels)  Variable 

Telescope  Primary  D  1.6  (m)  Constant 

Diameter 

Telescope  Secondary  D,  0.33  (m)  Constant 

Diameter 

Aperture  Radius  a  31  (pixels)  Constant 

Light  Wavelength  k  5,5  x  10'^  Constant 

(m) 

Cutoff  Frequency  f^  68.3  (1/arcsec)  Constant 

Random  Number  s  123456789  Constant 

Seed 
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arc-seconds 
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Figure  A3  Asteroid  Short  Exposure  Image: 

ro  =  0.206,  Np  =  10^. 


arc— seconds 

Figure  A.4  Asteroid  Short  Exposure  Image, 

ro  =  0.206,  Np  =  10*. 
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arc— seconds 

Figure  AS  Asteroid  Short  Exposure  Image, 


ro  =  0.206,  Np  =  10*. 


arc— seconds 

Figure  A.6  Asteroid  Short  Exposure  Image, 

ro  =  0.206,  Np  =  10*. 
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Figure  A.7  Asteroid  Short  Exposure  Image, 
To  =  0.103,  Np  =  10*. 


55 


arc— seconds  arc— seconds 


0.56 


-0.00 


-0.66  F 


-l.IO 


-1.10 


-0.55  -0.00  0.55 

arc-seconds 


Figiire  Asteroid  Short  Exposure  Image, 
ro  =  0.103,  N_  *  10*. 
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Figure  A.10  Asteroid  Short  Exposure  Image, 
ro  =  0.103,  Np  =  lO*. 
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Figure  A.11  Binary  Star  Short  Exposure  Image, 
ro  =  0.206,  Np  =  lO*. 


Figure  A.12  Binary  Star  Short  Exposure  Image, 

ro  =  0.206,  Np  =  10^, 
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Figure  A.18  Star  Short  Exposure  Image, 
ro  =  0.206,  Np  =  lO^. 
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Figure  A.19  Star  Short  Exposure  Image, 
ro  =  0.103,  Np  =  10*. 


Figure  A.20  Star  Short  Exposure  Image, 
ro  =  0.103,  Np  =  10^. 


61 


-1.10  -0.55  -0.00  0.55 

arc -seconds 

Figure  k22  Star  Short  Exposure  Image, 
ro  =  0.103,  Np  =  lO^. 
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Figure  A^<>  Baseline  KT  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  10*,  N,  =  400,  OS  =  4. 
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Figure  AJ7  Baseline  TC  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  10*,  N,  =  400,  OS  =  4. 
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RMS  P 


Figure  A^8  Asteroid  Baseline  Phasor  Spectrum  SNR. 


Spatial  Frequency  (l/arcsec) 

Figure  AJ29  Asteroid  Baseline  Azimuthal  RMS  Phase  Error. 
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Figure  A30  KT  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  10*,  N,  =  25,  OS  =  4. 


Figure  A31  TC  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  10*,  N,  =  25,  OS  =  4. 
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Figure  A32  Asteroid  Phaser  Spectrum  SNR,  Nf  =  25. 
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Figure  A33  Asteroid  Azimuthal  RMS  Phase  Error,  Nf  =  25. 
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Figure  A34  KT  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  10*,  N,  =  100,  OS  =  4. 
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Figure  AJ5  TC  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  10*,  Nj  =  100,  OS  =  4. 


69 


RMS  Phase  Error  (radians)  Phaser  Spectrum  SNR 


Figure  AJ6  Asteroid  Phasor  Spectrum  SNR,  Nf  =  100. 


Figure  A37  Asteroid  Azimuthal  RMS  Phase  Error,  Nf  =  100. 
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arc— seconds 

Figure  A38  KT  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  10*,  N,  =  1600,  OS  =  4. 
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Figure  AJ9  TC  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  10®,  N,  =  1600,  OS  =  4. 
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RMS  Phase  Error  (radians)  Phaser  Spectrum  SNR 


Figure  A.40  Asteroid  Phasor  Spectrum  SNR,  Nf  =  1600. 


Figure  A.41  Asteroid  Azimuthal  RMS  Phase  Error,  Nf  =  1600. 
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Figure  A.42  KT  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  10®,  N,  =  400,  OS  =  4. 
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gun  A.43  TC  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  10^,  N,  =  400,  OS  =  4. 
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RMS  Phase  Error  (radians)  Phaser  Spectrum  SNR 


Figure  A.44  Asteroid  Phaser  Spectrum  SNR,  Np  =  10*^. 


Figure  A.45  Asteroid  Azimuthal  RMS  Phase  Error,  Np  =  10^. 
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Figure  A.46  KT  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  10^,  N,  =  400,  OS  =  4. 
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Figure  A.47  TC  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  10^,  Nf  =  400,  OS  =  4. 
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RMS  Phase  Error  (radians)  Phaser  Spectrum  SNR 


Figure  A.48  Asteroid  Phaser  Spectrum  SNR,  Np  =  10^. 


Figure  A.49  Asteroid  Azimuthal  RMS  Phase  Error,  Np  =  10^. 
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Figure  AJO  KT  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  lO*,  N,  =  400,  OS  =  4. 


Figure  A.51  TC  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  lO’,  N^  =  400,  OS  =  4. 
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RMS  Phase  Error  (radians)  Phaser  Spectrum  SNR 


Figure  AJ2  Asteroid  Phaser  Spectrum  SNR,  Np  =  1(P. 


Figure  A.53  Asteroid  Azimuthal  RMS  Phase  Error,  Np  =  10^. 
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Figure  A.54  KT  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  10*,  N,  =  400,  OS  =  2. 
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Figure  A.55  TC  Recovered  Asteroid,  No  Tilt, 
ro  =  0.206,  Np  =  10*,  N,  =  400,  OS  =  2. 
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Figure  A.56  Asteroid  Phaser  Spectrum  SNR,  OS  =  2. 


Figure  A.57  Asteroid  Azimuthal  RMS  Phase  Error,  OS  =  2. 
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Figure  A^8  KT  Recovered  Asteroid,  No  Tilt, 
ro  =  0.103,  Np  =  10®,  Nf  =  400,  OS  =  2. 
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Figure  A.59  TC  Recovered  Asteroid,  No  Tilt, 
ro  =  0.103,  Np  =  10®,  Nf  =  400,  OS  =  2. 
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Figure  A.60  Asteroid  Phaser  Spectrum  SNR, 
ro  =  0.013,  OS  =  2. 


Spatial  Frequency  (1/arcsec) 

Figure  A.61  Asteroid  Azimuthal  RMS  Phase  Error, 
ro  =  0.103,  OS  =  Z 
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Figure  A.62  KT  Recovered  Asteroid, 
With  Tilt,  To  =  0.206,  N  =  10*, 

N,  =  400,  OS  =  4. 
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Figure  A.63  TC  Recovered  Asteroid, 
With  Tilt,  To  =  0.206,  N  =  10*, 

Ni  =  400,  OS  =  4. 
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Figure  kM  Asteroid  Phasor  Spectrum  SNR,  Tilt. 
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Figure  A.65  Asteroid  Azimuthal  RMS  Phase  Error,  Tilt 
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RMS  P 


Figure  A.68  Asteroid  Phaser  Spectrum  SNR, 
Tilt,  Np  *  10^. 


Spatial  Frequency  (1/arcsec) 

Figure  A.69  Asteroid  Azimuthal  RMS  Phase  Error, 
Tilt,  Np  *  10^. 
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arc- seconds  arc -seconds 


Figure  A.70  KT  Recovered  Asteroid, 
With  Tilt,  ro  =  0.206,  N  =  lO*, 

Nf  =  400,  OS  =  4. 


Figure  A.71  TC  Recovered  Asteroid, 
With  Tilt,  ro  =  0.206,  Np  =  1(>*, 

Nf  =  400,  OS  =  4. 
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Figure  A.72  Asteroid  Phaser  Spectrum  SNR, 

Tut,  Np  =  icy*. 


Figure  A.73  Asteroid  Azimuthal  RMS  Phase  Error, 
Tilt,  Np  =  10^. 
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Figure  A.74  KT  Recovered  Asteroid, 
With  Tilt,  ro  =  0.206,  N  =  10*, 

Nf  =  1600,  OS  =  4. 
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Figure  A.75  TC  Recovered  Asteroid, 
With  Tilt,  ro  =  0.206,  N  =  10*, 

Nt  =  1600,  OS  =  4. 
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Figure  A.78  KT  Recovered  Asteroid, 
With  Tilt,  Centroided,  ro  =  0.206, 
Np  =  10*,  N,  =  1600,  OS  =  4. 
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Figure  A.79  TC  Recovered  Asteroid, 
With  Tilt,  Centroided,  r©  =  0.206, 
Np  =  10*,  N,  =  1600,  OS  =  4. 
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Figure  A.80  Asteroid  Phaser  Spectrum  SNR, 
Tilt,  Centroiding,  Nf  =  1600. 


Figure  A.81  Asteroid  Azimuthal  RMS  Phase  Error,  Tilt, 
Centroiding,  Nf  =  1600. 
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Figure  A.82  KT  Recovered  Asteroid, 
With  Tilt,  Centroided,  r®  =  0.206, 
Np  =  10*,  Nf  =  400,  OS  =  4. 
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Figure  A.83  TC  Recovered  Asteroid, 
With  Tilt,  Centroided,  ro  =  0.206, 

Np  =  10*,  Nf  =  400,  OS  =  4. 
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RMS  Phase  Error  (radians)  Phaser  Spectrum  SNR 


Figure  A.84  Asteroid  Phaser  Spectrum  SNR,  Tilt,  Centroiding. 


Figure  A.85  Asteroid  Azimuthal  RMS  Phase  Error,  Tilt,  Centroiding. 
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Figure  A.86  KT  Recovered  Asteroid,  With  Tilt, 
Centroided,  r®  =  0.206,  N.  =  10^, 

Nf  =  400,  OS  =  4. 


arc-seconds 

Figure  A.87  TC  Recovered  Asteroid,  With  Tilt, 
Centroided,  ro  =  0.206,  N,  =  10*, 

Nf  =  400,  OS  =  4. 
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Phasor  Spectrum  SNR 


Figure  A-90  KT  Recovered  Asteroid,  With  Tilt, 
Centroided,  r^,  =  0.M6,  N.  =  lO*, 

N,  =  400,  OS  =  4. 


Figure  A.91  TC  Recovered  Asteroid,  With  Tilt, 
Centroided,  ro  =  0.206,  N-  =  10^, 

Nj  =  400,  OS  =  4. 
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RMS  Phase  Error  (radians)  Phaser  Spectrum  SNR 


Figure  A.96  Binaiy  Star  Phaser  Spectrum  SNR,  Np  =  1(P. 


Figure  A.97  Binaiy  Star  Azimuthal  RMS  Phase  Error,  Np  =  10^. 
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Figure  A.98  KT  Recovered  Binary  Star,  No  Tilt, 
ro  =  0.206,  Np  =  lO^,  N,  =  400,  OS  =  4. 


Figure  A.99  TC  Recovered  Binary  Star,  No  Tilt, 
ro  =  0.206,  Np  =  lO^,  Nf  =  400,  OS  =  4. 
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RMS  Phase  Error  (radians)  Phasor  Spectrum  SNR 


Figure  A.100  Binaiy  Star  Phasor  Spectrum  SNR,  Np  =  10^. 


Figure  A.101  Binary  Star  Azimuthal  RMS  Phase  Error,  Np  =  10^. 
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Figure  A.102  KT  Recovered  Binaiy  Star,  No  Tilt, 
To  =  0.103,  Np  =  IO3,  Nf  =  400,  OS  =  4. 


Figure  A.103  TC  Recovered  Binaiy  Star,  No  Tilt, 
ro  =  0.103,  Np  =  IO3,  N,  =  400,  OS  =  4. 
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Figure  A.106  KT  Recovered  Binary  Star,  No  Tilt, 
ro  =  0.103,  Np  =  10^,  Nf  =  400,  OS  =  4. 
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Figure  A.107  TC  Recovered  Binary  Star,  No  Tflt, 
ro  =  0.103,  Np  =  lO^,  N^  =  400,  OS  =  4. 
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Figure  A.108  Binary  Star  Phaser  Spectrum  SNR,  r^  =  0.103,  Np  =  lO' 


Figure  A.109  Binary  Star  Azimuthal  RMS  Phase  Error, 


APPENDK  B.  KNOX-THOMPSON  MAIN  PROGRAM 


c  THIS  PROGRAM  CREATES  ONE  OF  THREE  IMAGES:  A  STAR,  A 
c  BINARY  STAR  AND  AN  ASTEROID.  IT  THEN  DEGRADES  THE 
c  IMAGE  BY  SIMULATING  ATMOSPHERIC  CONDITIONS  AND  PHOTON 

c  NOISE,  AND  THEN  RECONSTRUCTS  THE  IMAGE  USING  THE 
c  KNOX-THOMPSON  ALGORITHM. 


c 

c 

c 

c 

c 

c 

c 

c 


AUTHOR:  LT  JAMES  M.  LACKEMACHER 

COMPL  DATE:  26  OCTOBER  1990 

REASON:  COMPLETE  REQUIREMENTS  FOR  A  MASTERS 

DEGREE  IN  PHYSICS. 

GOAL:  SIMULATE  OBJECT,  DEGRADE  OBJECT, 

RECONSTRUCT  OBJECT  USING  KNOX- 
THOMPSON  AND  TRIPLE-CORRELATION 
METHODS,  FILTER  AND  COMPARE. 


PROGRAM  KNOXTHOMPSON 


c  MAIN  PROGRAM  COMPLEX  VARIABLE  LIST 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


F  n  DIM  ARRAY  USED  IN  THE  FOURIER  TRANSFORM 

GAUSSIAN  n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  GAUSSIAN 
PORTION  OF  THE  ASTEROID 
I  n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE 

DEGRADED  IMAGE 

nCT  nxnx5x9  ARRAY  THAT  REPRESENTS  THE  CROSS¬ 

SPECTRUM  OF  THE  IMAGE 

IS  n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE 

DEGRADED  POINT  SOURCE 

ISKT  n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  MODULUS 

SQUARED  OF  THE  POINT  SOURCE 
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c 

c 

c 

c 

c 

c 

c 


c 


OKT 

OBJDATA 

KTPHASOR 

PUPIL 

TEMPDATA 


n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  OBJECT 
SPECTRUM  OF  THE  RECONSTRUCTED  IMAGE 
n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  OBJECT 
n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  PHASOR  OF 
THE  OBJECT  IN  THE  RECONSTRUCHON  PROCESS 
n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  PUPIL 
PORTION  OF  THE  ASTEROID 
n  X  n  DIM  ARRAY  THAT  IS  USED  AS  A  TEMPORARY 
ARRAY 


c 


MAIN  PROGRAM  REAL  VARIABLE  UST 


c  KTsnr 

c 

c  mod 
c 

c  rsnr 
c 

c  xvarrkt 

c 

c 

c  xvarikt 

c 

c 


n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  SNR  OF 
EACH  PHASOR 

n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  MODULUS 
OF  THE  RECONSTRUCTED  IMAGE 
n/2  DIM  ARRAY  THAT  REPRESENTS  THE  SNR  AS  A 
FUNCTION  OF  RADIUS 

n  X  n  X  5  X  9  DIM  ARRAY  THAT  REPRESENTS  THE 
REAL  PART  OF  THE  VARIANCE  OF  THE  CROSS- 
SPECTRUM 

n  X  n  X  5  X  9  DIM  ARRAY  THAT  REPRESENTS  THE 
IMAGINARY  PART  OF  THE  VARIANCE  OF  THE  CROSS- 
SPECTRUM 


c  MAIN  PROGRAM  INTEGER  VARIABLE  UST 


c  fwd 
c  icounter 
c 

c  inv 
c  In 
c  o^et 
c 
c 

c  mseed 
c 


VALUE  OF  1  FOR  FORWARD  FFT 
COUNTER  THAT  COUNTS  THE  NUMBER  OF 
SNAPSHOTS 

VALUE  OF  -1  FOR  INVERSE  FFT 
2^^^  FOR  USE  WITH  FFT  SUBROUTINE 
VARIABLE  THAT  REPRESENTS  THE  NUMBER  OF 
PIXELS  THAT  ARE  AVERAGED  IN  THE  KNOX- 
THOMPSON  PROCESS 

VARIABLE  USED  TO  ENSURE  ONLY  ONE  PASS  OF 
INITIAL  PART  OF  PHSUB  SUBROUTINE 
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c  n 
c  nframes 
c  nphoton 
c  nyquist 
c 
c 
c 


DIMENSION  OF  ONE  SIDE  OF  2-DIM  ARRAY 
TOTAL  NUMBER  OF  SHORT  EXPOSURE  SNAPSHOTS 
TOTAL  NUMBER  OF  PHOTONS  IN  SNAPSHOT 
EQUAL  TO  THE  TELESCOPE  PUPIL  FUNCTION  RADIUS 
DERIVED  FROM  THE  RELATION: 
nyquist  =  (telescope  diameter  x 
number  of  pixels  per  rO)/rO 


MAIN  PROGRAM 


PARAMETER(ii=64,ln=6,fwd-  l,inv=-l) 

COMPLEX*  16  OBJDATA(n,n),  TEMPDATA(n,n),  F(n), 

+  PUPIL(n,n),  STARDATA(n,n),  GAUSSIAN(n,n), 

+  KTPHASOR(n,n),  IKT(n,n,5,9),  ISKT(n,n), 

+  I(n,n),  IS(n,n),  OKT(n,n) 

REAL*8  mod(n,n),  xvarrkt(n,n,5,9),  xvarilct(n,n,5,9), 

+  KTsnr(n,n),  rsnr(n/2) 

INTEGER  offset 

CHARACTER*  16  filel,  file2 
CHARACTER*!  cent 

c  INITIALIZE  MSEED  TO  ALLOW  ONLY  ONE  PASS  THROUGH  FIRST 
c  PART  OF  PHSUB 

mseed  =  1 

c  INITIALIZE  PROGRAM  READING  REQUIRED  VARIABLES  AND 
c  CREATE  THE  DESIRED  OBJECT 

CALLInitialize(OBJDATA,GAUSSIAN,PUPIL,F,TEMPDATA, 

+  filel, file2,fwd,inv, In, nframes, nphoton, nyquist, 

+  offset,cent,n) 

c  TRANSFORM  THE  OBJECT  TO  FREQUENCY  SPACE 
CALL  FFT2D(OBJDATA,TEMPDATA,F,ln,fwd,n) 
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c  NOR^MLIZE  ALL  OF  THE  ARRAY  ELEMENTS  TO  THE  VALUE  OF 
c  THE  NUMBER  OF  PHOTONS  WHERE  THE  DC  TERM  EQUALS  THE 
c  NUMBER  OF  PHOTONS  OF  THE  IMAGE 

CALL  Photons(OBJDATA,nphoton^) 

c  CREATE  THE  POINT  SOURCE 

CALL  Star(STARDATA,n) 

c  TRANSFORM  THE  POINT  SOURCE  TO  FREQUENCY  SPACE 

CALL  FFT2D(STARDATA,TEMPDATA,F4n,fwd,n) 

c  NORMALIZE  ALL  OF  THE  ARRAY  ELEMENTS  TO  THE  VALUE  OF 
c  THE  NUMBER  OF  PHOTONS  WHERE  THE  DC  TERM  EQUALS  THE 
c  NUMBER  OF  PHOTONS  OF  THE  IMAGE 

CALL  Photons(STARDATA,nphoton,n) 

c  COMMENCE  THE  LOOP  THAT  COUNTS  THE  SNAPSHOTS 

DO  10  icounter  «  1,  nfirames 

c  DEGRADE  THE  IMAGE  WITH  THE  TELESCOPE,  THE  ATMOSPHERE, 
c  AND  THE  PHOTON  NOISE 

CALL  PHSUB(I,F,TEMPDATA,OBJDATA,icounter,nframes, 

+  nyquist,fwd,inv,ln,n,inseed) 

c  CENTROID  THE  DEGRADED  IMAGE  ONLY  SINCE  ONLY  PHASE  IS 
c  EFFECTED  BY  CENTROIDING  IF  CENTROIDING  IS  DESIRED. 

IF  ((ccnt.EQ.’Y’).OR.(cent.EQ.y))  THEN 
CALL  Centroid(I,TEMPDATA,n) 

ENDIF 

c  TRANSFORM  THE  DEGRADED  IMAGE  TO  FREQUENCY  SPACE 
CALL  FFT2D(I,TEMPDATA,F,ln,fwd,n) 


no 


c  DEGRADE  THE  POINT  SOURCE  WITH  THE  TELESCOPE,  THE 
c  ATMOSPHERE,  AND  THE  PHOTON  NOISE 

CALLPHSUB(IS,F,TEMPDATA,STARDATA,icounter,nframes, 

+  nyquist,fwd,inv,ln,n,inseed) 

c  TRANSFORM  THE  DEGRADED  IMAGE  TO  FREQUENCY  SPACE 

CALL  FFT2D(IS,TEMPDATA,F,ln,fwd,n) 

c  DETERMINE  THE  MODULUS  OF  THE  DEGRADED  IMAGE  AND 
c  POINT  SOURCE 

CALL  Modulus(I,IS,IKT,ISKT,mocl,icounter, 

+  nfirames,nphoton,n) 

c  RECONSTRUCT  THE  OBJECT  FROM  THE  DEGRADED  IMAGE  AND 
c  POINT  SOURCE 

CALL  KTrecon(KTPHASOR,I,IKT,KTsnr,offset,xvarrkt, 

+  xvarikt,icounter,nfranies,nyquist,n) 

WRITE(*,*)icounter,TRAMES  COMPLETED’ 

c  END  THE  LOOP 

10  CONTINUE 

c  CALCULATE  THE  AVERAGE  SNR  AS  A  FUNCTION  OF  RADIUS 

CALL  SNRcalc(KTsnr,rsnr,nyquist,n) 

c  COMBINE  THE  MODULUS  WITH  THE  PHASOR 

CALL  Combine(OKT,raod,KTPHASOR,n) 

c  WRITE  RECONSTRUCTED  PHASE  AND  POWER  SPECTRUM  TO  A 
c  FILE 

CALL  Writefile(OKT,rsnr,file  l,file2,nyquist,n) 

STOP 

END 
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APPENDIX  C  TRIPLE-CORRELATION  MAIN  PROGRAM 


c  THIS  PROGRAM  CREATES  ONE  OF  THREE  IMAGES:  A  STAR,  A 
c  BINARY  STAR  AND  AN  ASTEROID.  IT  THEN  DEGRADES  THE 

c  IMAGE  BY  SIMULATING  ATMOSPHERIC  CONDITIONS  AND  PHOTON 

c  NOISE,  AND  THEN  RECONSTRUCTS  THE  IMAGE  USING  THE  TRIPLE- 

c  CORRELATION  ALGORITHM. 


c 

c 

c 

c 

c 

c 

c 

c 


AUTHOR:  LT  JAMES  M.  LACKEMACHER 

COMPL  DATE:  26  OCTOBER  1990 

REASON:  COMPLETE  REQUIREMENTS  FOR  A  MASTERS 

DEGREE  IN  PHYSICS. 

GOAL:  SIMULATE  OBJECT,  DEGRADE  OBJECT, 

RECONSTRUCT  OBJECT  USING  KNOX- 
THOMPSON  AND  TRIPLE-CORRELATION 
METHODS,  FILTER  AND  COMPARE. 


PROGRAM  TRIPLECORR 


c  MAIN  PROGRAM  COMPLEX  VARIABLE  LIST 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


BSPHASOR 

F 

GAUSSIAN 

I 

IBS 

IDBS 


n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  PHASOR  OF 
THE  OBJECT  IN  THE  RECONSTRUCTION  PROCESS 
n  DIM  ARRAY  USED  IN  THE  FOURIER  TRANSFORM 
n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  GAUSSIAN 
PORTION  OF  THE  ASTEROID 

n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  DEGRADED 
IMAGE 

n  X  n  X  5  X  9  ARRAY  THAT  REPRESENTS  THE 
BISPECTRUM  OF  THE  IMAGE 
5x9  ARRAY  THAT  REPRESENTS  THE  DEVIATION 
FROM  THE  DC  TERM 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


IS 

ISBS 

OBS 

OBJDATA 

PUPIL 

TEMPDATA 


n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  DEGRADED 
POINT  SOURCE 

n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  MODULUS 
SQUARED  OF  THE  POINT  SOURCE 
n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  OBJECT 
SPECTRUM  OF  THE  RECONSTRUCTED  IMAGE 
n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  OBJECT 
n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  PUPIL 
PORTION  OF  THE  ASTEROID 
n  X  n  DIM  ARRAY  THAT  IS  USED  AS  A  TEMPORARY 
ARRAY  IN  THE  FOURIER  TRANSFORM 


c 


MAIN  PROGRAM  REAL  VARIABLE  LIST 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


BSsnr 

mod 

rsnr 

xvarrbs 

xvaribs 


n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  SNR  OF 
EACH  PHASOR 

n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE 
MODULUS  OF  THE  RECONSTRUCTED  IMAGE 
n/2  DIM  ARRAY  THAT  REPRESENTS  THE  SNR  AS  A 
FUNCTION  OF  RADRJS 

n  X  n  X  5  X  9  DIM  ARRAY  THAT  REPRESENTS  THE 
REAL  PART  OF  THE  VARIANCE  OF  THE  BISPECTRUM 
n  X  n  X  5  X  9  DIM  ARRAY  THAT  REPRESENTS  THE 
IMAGINARY  PART  OF  THE  VARIANCE  OF  THE 
BISPECTRUM 


c  MAIN  PROGRAM  INTEGER  VARIABLE  LIST 


c  fwd 

c  icounter 

c 

c  inv 

c  In 

c  offset 


VALUE  OF  1  FOR  FORWARD  FFT 
COUNTER  THAT  COUNTS  THE  NUMBER  OF 
SNAPSHOTS 

VALUE  OF  -1  FOR  INVERSE  FFT 
2 "'In  FOR  USE  WITH  FFT  SUBROUTINE 
VARIABLE  THAT  REPRESENTS  THE  NUMBER  OF 
PIXELS  THAT  ARE  AVERAGED  IN  TFE  TRIPLE¬ 
CORRELATION  PROCESS 
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+  + 


c  mseed 
c 

c  n 
c  nframes 
c  nphoton 
c  nyquist 
c 


VARIABLE  USED  TO  ENSURE  ONLY  ONE  PASS  OF 
INITIAL  PART  OF  PHSUB  SUBROUTINE 
DIMENSION  OF  ONE  SIDE  OF  2-DIM  ARRAY 
TOTAL  NUMBER  OF  SHORT  EXPOSURE  SNAPSHOTS 
TOTAL  NUMBER  OF  PHOTONS  IN  SNAPSHOT 
EQUAL  TO  THE  TELESCOPEPUPIL  FUNCTION  RADIUS 
DERIVED  FROM  THE  RELATION: 


c  nyquist  =  (telescope  diameter  x 

c  number  of  pixels  per  rO)/iO 


MAIN  PROGRAM 


PARAMETER(n=64,ln=6,fwd= l,inv=-l) 


COMPLEX*16  OBJDATA(n,n),  TEMPDATA(n,n),  F(n), 

+  PUPIL(n,n),  I(n,n),  STARDATA(n,n),  GAUSSIAN(n,n), 

+  BSPHASOR(n,n),  IBS(n4i,5,9),  ISBS(n,n),  IS(n,n), 

+  OBS(n,n),  IDBS(5,9) 

REAL*8  mod(n,n),  xvarrbs(n^,5,9),  xvaribs(n,n,5,9), 

+  BSsnr(n,n),  rsnr(n/2) 

INTEGER  offset 

CHARACTER*16  fflel,  file2 
CHARACTER*  1  cent 

c  DSrmALIZE  MSEED  TO  ALLOW  ONLY  ONE  PASS  THROUGH  FIRST 
c  PART  OF  PHSUB 

mseed  =  1 

c  INITIALIZE  PROGRAM  READING  REQUIRED  VARIABLES  AND 
c  CREATE  THE  DESIRED  OBJECT 

CALLInitialize(OBJDAT>^GAUSSIAN,PUPIL,F,TEMPDATA, 
filel,file2,fwd,inv,ln,inframes, 
nphoton,nyquist,offseLcent,n) 
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c  TRANSFORM  THE  OBJECT  TO  FREQUENCY  SPACE 

CALL  FFT2D(OBJDATA,TEMPDATA,F,ln4wd,n) 

c  NORMALIZE  ALL  OF  THE  ARRAY  ELEMENTS  TO  THE  VALUE  OF 
c  THE  NUMBER  OF  PHOTONS  WHERE  THE  DC  TERM  EQUALS  THE 
c  NUMBER  OF  PHOTONS  OF  THE  IMAGE 

CALL  Photons(OBJDATA,nphotoiMi) 

c  CREATE  THE  POINT  SOURCE 

CALL  Star(STARDATA,n) 

c  TRANSFORM  THE  POINT  SOURCE  TO  FREQUENCY  SPACE 

CALL  FFT2D(STARDATA,TEMPDATA,F4n,fwd,n) 

c  NORMALIZE  ALL  OF  THE  ARRAY  ELEMENTS  TO  THE  VALUE  OF 
c  THE  NUMBER  OF  PHOTONS  WHERE  THE  DC  TERM  EQUALS  THE 
c  NUMBER  OF  PHOTONS  OF  THE  IMAGE 

CALL  Photons(STARDATA,nphoton,n) 

c  COMMENCE  THE  LOOP  THAT  COUNTS  THE  SNAPSHOTS 

DO  10  icounter  =  1,  nframes 

c  DEGRADE  THE  IMAGE  WITH  THE  TELESCOPE,  THE  ATMOSPHERE, 

c  AND  THE  PHOTON  NOISE 

CALL  PHSUB(I,F,TEMPDATA,OBJDATA,icounter,nfirames, 

+  nyquist,fwd,inv,ln,n,mseed) 

c  CENTROID  THE  DEGRADED  IMAGE  ONLY  SINCE  ONLY  PHASE  IS 
c  EFFECTED  BY  CENTROIDING  IF  CENTROIDING  IS  DESIRED. 

IF  ((cent.EQ.’Y’).OR.(cent.EQ.y))  THEN 
CALL  Centroid(I,TEMPDATA,n) 

ENDIF 
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+  + 


c  TRANSFORM  THE  DEGRADED  IMAGE  TO  FREQUENCY  SPACE 

CALL  FFT2D(I,TEMPDATA,F.ln4wd,n) 

c  DEGRADE  THE  POINT  SOURCE  WITH  THE  TELESCOPE,  THE 
c  ATMOSPHERE,  AND  THE  PHOTON  NOISE 

CALL  PHSUB(IS,F,TEMPDATA,STARDATA,icounter,nframes, 

+  nyquist,fwd,inv,ln,n,mseed) 

c  TRANSFORM  THE  DEGRADED  IMAGE  TO  FREQUENCY  SPACE 

CALL  FFT2D(IS,TEMPDATA,F4n4wd,n) 

c  DETERMINE  THE  MODULUS  OF  THE  DEGRADED  IMAGE  AND 

c  POINT  SOURCE 

CALL  Modulus(I,IS,IBS,ISBS4nod,icounter, 

+  nframes,nphoton,n) 

c  RECONSTRUCT  THE  OBJECT  FROM  THE  DEGRADED  IMAGE  AND 
c  POINT  SOURCE 

CALL  BSrecon(BSPHASOR,I,IBS,IDBS,BSsnr,offset, 
xvarrbs,xvaribs,icounter,n£raines, 
nyquist,nphoton,n) 

WRITE(*,*)icounter,TRAMES  COMPLETED’ 
c  END  THE  LOOP 
10  CONTINUE 

c  CALCULATE  THE  AVERAGE  SNR  AS  A  FUNCTION  OF  RADIUS 

CALL  SNRcalc(BSsnr,rsnr,nyquist,n) 
c  COMBINE  THE  MODULUS  WITH  THE  PHASOR 
CALL  Combine(OBS,mod,BSPHASOR,n) 
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c  WRITE  RECONSTRUCTED  PHASE  AND  POWER  SPECTRUM  TO  A 

c  FILE 


CALL  Writefile(OBS^nr,filel^e2,]:7quist4i) 

STOP 

END 
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APPENDK  D.  UNIVERSAL  IMAGE  RECOVERY  SUBROUTINES 


c  THE  FOLLOWING  SUBROUTINES  WERE  GENERATED  BY  THE 

c  THESIS  AUTHOR  AND  ARE  REQUIRED  BY  BOTH  THE  KNOX- 
c  THOMPSON  AND  THE  TRIPLE-CORRELATION  PROGRAMS, 
c  ADDITIONALLY,  SOME  SUBROUTINES  ARE  REQUIRED  FOR 
c  THE  FILTRMS  PROGRAM  IN  APPENDIX  H. 


SUBROUTINE  LIST 


SUBROUTINE  Ast(ASTEROID, GAUSSIAN, PUPIL,F,TEMPDATA, 
+  fwd,inv4n,n) 

COMPLEX*16  GAUSSIAN(n4i),  ASTEROID(n,n),  PUPIL(n,n), 

+  TEMPDATA(n,n),  F(ii),  El,  E2,  E3,  E4,  E5,  E6,  E7 
maxval  »  0.0 
aval  =  4.0 
n2pl  =  11/2+1 
rad  =  16.0 
DO  10  i  «  1,  n 
DO  10  j  =  l,n 
X  =  float(j  -  (n2pl)) 
y  =  £loat((n2pl)  -  i) 
radius  =  sqrt(x**2.0  +  y**2.0) 

IF  (radius.LE.rad)  THEN 
PUPIL(y)  =  (1.0,0.0) 

ELSE 

PUPIL(id)  =  (0.0,0.0) 

ENDIF 

IF  ((abs(x).LE.aval).AND.(abs(y).LE.aval))  THEN 
GAUSSIAN(id)  =  DCMPLX(exp(-(x**2.0  + 

+  y**2.0)/4)) 

ELSE 

GAUSSIAN(iJ)  =  (0.0,0.0) 

ENDIF 

10  CONTINUE 

CALL  FFT2D(GAUSSIAN,TEMPDATA,F,ln,fwd,n) 

CALL  FFT2D(PUPIL,TEMPDATA,F,In,fwd,n) 
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DO  20  i  =  1,  n 
DO20j  =  l,n 

ASTEROID(ij)  =  GAUSSIAN(io)  *  PUPIL(ij) 

20  CX)NTINIJE 

CALL  FFT2D(ASTEROID,TEMPDATAF,ln,mv,n) 
a  =  10 
al  =  2.5 
a2  =  2 
a3  =  2.5 
a4  =  2 
a5  =  3.5 
a6  =  3 
a?  =  2.5 
DO  30  i  =  1,  n 
DO  30  j  =  1,  n 
X  =  float(j  -  (n2pl)) 
y  =  float((n2pl)-  i) 
xl  =  x+6 
yl  =  y-8 
xal  =  xl/al 
yal  =  yl/al 

IF  ((abs(xal).LE.aval)AND.(abs(yal).LE.aval)) 

+  THEN 

El  =  DCMPLX(a  *  exp(.(xal**2.0  +  yal**2.0))) 
ELSE 

El  =  (0.0, 0.0) 

ENDIF 
x2  =  x-4 
y2  =  y-6 
xa2  =  x2/a2 
ya2  =  y2/a2 

IF  ((abs(xa2).LE.aval)AND.(abs(ya2).LE.aval)) 

+  THEN 

E2  =  DCMPLX(a  *  exp(-(xa2**2.0  +  ya2**2.0))) 
ELSE 

E2  =  (0.0, 0.0) 

ENDIF 
x3  =  x-10 
y3  =  y-4 
xa3  =  x3/a3 
ya3  =  y3/a3 
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IF  ((abs(xa3).LE.aval)AND.(abs(ya3).LE.aval)) 

+  THEN 

E3  =  DCMPLX(a  *  exp(-(xa3**2.0  +  ya3**2.0))) 
ELSE 

E3  =  (0.0, 0.0) 

ENDIF 
x4  =  x+6 
y4  =  y+6 
xa4  «  x4/a4 
ya4  «  y4/a4 

IF  ((abs(xa4).LE.avaI).AND.(abs(ya4).L£.aval)) 

+  THEN 

E4  =  DCMPLX(a  *  exp(-(xa4**2.0  +  ya4**2.0))) 
ELSE 

E4  =  (0.0, 0.0) 

ENDIF 
x5  =  x+0 
y5  =  y.2 
xa5  =  x5/a5 
ya5  =  y5/a5 

IF  ((abs(xa5).LE.aval).AND.(abs(ya5).LE.aval)) 

+  THEN 

E5  =  DCMPLX(a  *  exp(-(xa5**2.0  +  ya5**2.0))) 
ELSE 

E5  =  (0.0,0.0) 

ENDIF 
x6  =  x+2 
y6  =  y+8 
xa6  =  x6/a6 
ya6  =  y6/a6 

IF  ((abs(xa6).L£.aval).AND.(abs(ya6).LE.ava])) 

+  THEN 

E6  =  DCMPLX(a  •  exp(-(xa6**2.0  +  ya6**2.0))) 
ELSE 

E6  =  (0.0,0.0) 

ENDIF 
x7  =  x-6 
y7  =  y+6 
xa7  =  x7/a7 
ya7  =  y7/a7 
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IF  ((abs(xa7).LE.aval)j\ND.(abs(ya7).LE.aval)) 

+  THEN 

E7  =  DCMPLX(a  •  exp(-(xa7**2.0  +  ya7**2.0))) 
ELSE 

E7  =  (0.0, 0.0) 

ENDIF 

ASTEROID(ij)  =  DCMPLX(DREAL(ASTEROID(y)  - 
+  ^1  +  E2  +  E3  +  E4  +  E5  +  E6  +  E7))) 

IF  (DREAL(ASTEROID(ij)).LT.O.O)  ASTEROID(ij)  = 

+  (0.0, 0.0) 

30  CONTINUE 
RETURN 
END 


SUBROUTINE  Bistar(DATA,n) 

c  THIS  S/R  CREATES  A  SIMULATED  BINARY  STAR  WITH  ONE  STAR 
c  LARGER  THAN  THE  OTHER 

COMPLEX*16  DATA(n,n) 
n2pl  =  n/2  +  1 
DO  10  i  =  1,  n 
IX)10j  =  l,n 
X  *  float(j  -  n2pl) 
y  =  float(n2pl  -  i) 

IF  ((x.EQ.12.0).AND.(y.EQ.12.0))  THEN 
DATA(ij)  =  (2.0,0.0) 

ELSEIF  ((x.EQ.-12.0)j\ND.(y.EQ.-12.0))  THEN 
DATA(ij)  =  (1.0,0.0) 

ELSE 

DATA(ij)  =  (0.0, 0.0) 

ENDIF 

10  CONTINUE 
RETURN 
END 


121 


SUBROUTINE  Centroid(DATA,TEMPDATA,n) 

c  THIS  S/R  DETERMINES  THE  CENTROID  OF  THE  DEGRADED  IMAGE 

c  AND  THEN  CENTROIDS  THE  IMAGE 

COMPLEX*16  DATA(n,n),  TEMPDATA(n,n) 

REAL*8  yi,  xnum,  ynum,  xden,  yden 
n2pl  =  11/2+1 
DO  10  i  =  1,  n 
DO  10  j  =  1,  n 
ig  =  dfloatQ  -  n2pl) 
yi  =  d£Ioat(ii2pl  •  i) 
xnum  =  xnum  +  jg  *  ABS(DATA(ij)) 
ynum  =  ynum  +  yi  *  ABS(DATA(ij)) 
xden  =  xden  +  ABS(DATA(y)) 
yden  =  yden  +  ABS(DATA(ij)) 

10  CONTINUE 

jxbar  =  idnint(xnum^den) 
iybar  =  idnint(ynum^den) 

DO  20  i  =  1,  n 
DO20j  =  l,n 
ii  =  i  -  iybar 
jj  =  j  +  jxbar 

IF((ii.GT.n).OR.(ii.LT.l).OR.(ij.GT.n).OR. 

+  (ij.LT.l))  THEN 

TEMPDATA(ij)  =  (0.0, 0.0) 

ELSE 

TEMPDATA(ij)  =  DATA(iijj) 

ENDIF 

20  CONTINUE 
DO  30  i  =  1,  n 
DO30j  =  l,n 

DATA(io)  =  TEMPDATA(y) 

30  CONTINUE 
RETURN 
END 
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SUBROUTINE  Combme(0,mod,PHASOR,n) 

c  THIS  S/R  CX)MBINES  THE  MODULUS  AND  PHASOR  OF  THE 
c  RECONSTRUCTED  OBJECT 

COMPLEX*16  PHASOR(-ny2:ny2-l,-iiy2:n/2-l), 

+  0(-n/2:n/2-l,-n/2:n/2-l) 

REAL*8  mod(-n/2:n/2-l,-n/2:n/2-l) 
n2  =  ii/2 
n2ml  =  ii2  - 1 
DO  10  i  =  -n2,  n2ml 
DO  10  j  =  -Til,  n2ml 
0(y)  =  mod(ij)  *  PHASOR(ij) 

10  CONTINUE 
RETURN 
END 

SUBROUTINE  Complexconj(DATA,n) 

c  THIS  S/R  IS  CALLED  BY  FFT2D  AND  DETERMINES  THE  COMPLEX 
c  CONJUGATE  OF  THE  2-D  ARRAY  IN  ORDER  FOR  THE  ARRAY  TO 
c  BE  INVERSE  FFTed. 

COMPLEX*16  DATA(n,n) 

DO  10  i  =  1,  n 
DO  10  j  =  1,  n 

DATA(id)  =  DCONJG(DATA(id)) 

10  CONTINUE 
RETURN 
END 
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SUBROUTINE  FFT(F,ln,n) 

c  THIS  S/R  IS  USED  BY  THE  2-D  FFT  TWICE,  FIRST  FFTing  THE 
c  ROWS  THEN  THE  COLUMNS  OF  THE  2-D  ARRAY.  THIS  S/R  WAS 
c  AQUIRED  FROM  'DIGITAL  IMAGE  PROCESSING"  BY  GONZALEZ 
c  AND  WINTZ. 

COMPLEX*16  F(n),U,W,T 
REAL*8  pi,  one 
one  =  l.OD+00 
pi  =  DACOS(-one) 
nv2  =  n/2 
nml  =  n  - 1 

j  =  1 

DO  3  i  =  1,  nml 
IF(i.GE,j)  GOTO  1 
T  =  Fa) 
m  =  F(i) 

F(i)  =  T 

1  k  =  nv2 

2  IF(k.GE.j)  GOTO  3 
j  =  j-k 

k  =  k/2 
GOTO  2 

3  j  =  j  +  k 
DO  5  1  =  1,  In 

le  =  2**1 
lei  =  le/2 
U  =  (1.0,0.0) 

W  =  DCMPLX(DCOS(pi/lel),-DSIN(pi/lel)) 

DO  5  j  =  1,  lei 
DO  4  i  =  j,n,le 
ip  =  i  +  lel 
T  =  F(ip)  *  U 
F(ip)  =  F(i)  -  T 

4  F(i)  =  F(i)  +  T 

5  U  =  U  •  W 
RETURN 
END 
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SUBROUTINE  FFT2D(DATA,TEMPDATA,F,ln,dir,n) 

c  THIS  S/R  IS  THE  MAIN  S/R  AND  PERFORMS  A  2-D  FFT  ON  AN 
c  ARRAY  OF  SIDE  LENGTH  ’n’  WHERE  Tn’  IS  THE  POWER  OF  2. 
c  ’dir’  IS  EITHER  + 1  OR  -1  WHETHER  XFORMING  OR  INVERSE 

c  XFORMING  RESPECTIVELY.  ’TEMPDATA’ IS  A  WORKING  ARRAY 
c  FOR  THE  QUADRANT  SWAPPING  S/R.  ’F  IS  THE  1-D  ARRAY  USED 
c  BY  THE  1-D  FFT  S/R.  THIS  S/R  CALLS  QUADSWAP,  FFT,  FOR  BOTH 
c  FORWARD  AND  INVERSE  FFT  AND  CALLS  NORMFFT  AND 
c  COMPLEXCONJ  FOR  INVERSE  FFTs  ONLY.  THIS  FFT  NORMALIZES 
c  BY  DIVIDING  BY  n**2  WHEN  THE  INVERSE  FFT  IS  PERFORMED. 

COMPLEX*16  DATA(n,n),  F(n),  TEMPDATA(n,n) 

INTEGER  dir 

CALL  Quadswap(DATA,TEMPDATA,n) 

IF  (dir.EQ.-l)  THEN 
CALL  Complexconj(DATA,n) 

ENDEF 

DO  10  i  =  1,  n 
DO20j  =  l,n 
F(j)  *  DATA(id) 

20  CONTINUE 
CALL  FFT(F,ln,n) 

DO  30  j  «  1,  n 
DATA(id)  =  FO) 

30  CONTINUE 
10  CONTINUE 
DO  40  j  =  1,  n 
DO  50  i  =  1,  n 
F(i)  =  DATA(ij) 

50  CONTINUE 
CALL  FFT(F,ln,n) 
do  60  i  =  1,  n 
DATA(id)  =  F(i) 

60  CONTINUE 
40  CONTINUE 

IF  (dir.EQ.-l)  THEN 
CALL  Complexconj(DATA,n) 

CALL  NormFFT(DATAn) 

ENDIF 

CALL  Quadswap(DATA,TEMPDATA,n) 

RETURN 

END 
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SUBROUTINE  Modulus(I,IS,IDATA,ISDATA,mod,icounter, 

+  nfiraines,nphoton,n) 

c  THIS  S/R  DETERMINES  THE  MODULUS  OF  THE  DEGRADED  IMAGE 
c  AND  THE  POINT  SOURCE  AND  NORMALIZES  THEM  BY  DIVIDING 
c  BY  THEIR  RESPECTIVE  DC  VALUES. 


COMPLEX*  16  IDATA(-n/2:n/2-l,-n/2:n/2-l,0:4,-4:4), 
+  ISDATA(-n/2:n/2-l,-n/2:n/2-l),  DC,  DCS, 

+  I(-n/2:n/2-l,.ii/2:n/2-l), 

+  IS(-n/2:n/2-l,-n/2:n/2-l) 

REAL*8  mod(-ii/2:n/2-l,-n/2:n/2-l) 
k  =  0 
ii2  =  n/2 
n2ml  =  n2  - 1 
DO  10  ii  =  -n2,  n2ml 
DO  10  jj  =  -n2,  n2ml 

IDATA(iijj,k,k)  =  IDATA(iijj,k,k)  +  ((I(iijj)  * 

+  DCONJG(I(iiji)))  -  I(k,k)) 

ISDATA(ujj)  =  ISDATA(iiaj)  +  ((IS(iijj)  * 

+  DCONJG(IS(iijj)))  -  IS(k,k)) 

10  CONTINUE 

IF  (icounter.EQ.nframes)  THEN 
DC  =  IDATA(k,k,k,k) 

DCS  =  ISDATA(k,k) 

DO  20  ii  =  -n2,  n2ml 
DO  20  jj  =  -n2,  n2ml 
IDATA(iijj,k,k)  =  IDATA(iijj,k,k)/DC 
ISDATA(iijj)  =  ISDATA(iiaj)/DCS 
niod(iijj)  =  nphoton  * 

+  dsqrt(ABS(IDATA(ujj,k,k)/ISDATA(iioj))) 
IF  (mod(ii  jj).GT.dfloat(nphoton)) 

+  mod(iijj)  =  dfloat(nphoton) 

20  CONTINUE 
ENDIF 
RETURN 
END 
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SUBROUTINE  NormFFT(DATA,n) 

c  THIS  S/R  IS  CALLED  BY  FTT2D  FOR  INVERSE  FOURIER  XFORMS 
c  AND  NORMALIZES  THE  XFORM  BY  DIVIDING  BY  n**2. 

COMPLEX*  16  DATA(n,n) 
nsqrd  =  n  *  n 
DO  10  i  =  1,  n 
DO  10  j  =  1,  n 
DATA(iJ)  =  DATA(ij)/nsqrd 
10  CONTINUE 
RETURN 
END 

SUBROUTINE  Photons(DATA,nphoton,n) 

c  THIS  S/R  NORMALIZES  THE  DATA  ARRAY  TO  MAKE  THE  DC 
c  VALUE  EQUAL  TO  THE  NUMBER  OF  PHOTONS  IN  THE  OBJECT. 

COMPLEX*  16  DATA(n,n) 

REAL*8  photonum 
ic  =  n/2  +  1 

photonum  =  dfloat(nphoton)/DREAL(DATA(ic,ic)) 

DO  10  i  *  1,  n 
DO  10  j  =  1,  n 

DATA(ij)  =  photonum  *  DATA(ij) 

10  CONTINUE 
RETURN 
END 


SUBROUTINE  Quadswap(DATA,TEMPDATA,n) 

c  THIS  S/R  IS  CALLED  BY  FFT2D  AND  SWAPS  QUADRANTS  OF  DATA 
c  ARRAY  USING  TEMPDATA  AS  A  WORKING  ARRAY. 

COMPLEX*  16  DATA(n,n),  TEMPDATA(n,n) 
n2  =  n/2 
DO  10  i  =  1,  n/2 
DO  10  j  =  1,  n/2 

TEMPDATA(ij)  =  DATA(i+n2j+n2) 

TEMPDATA(i+n2j)  =  DATA(i,j+n2) 

TEMPDATA(ij+n2)  =  DATA(i+n2j) 
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TEMPDATA(i+n2o+ii2)  =  DATA(ij) 
10  CONTINUE 
DO  20  i  =  1,  n 
DO20j  =  1,  n 

DATA(ij)  =  TEMPDATA(y) 

20  CONTINIJE 
RETURN 
END 


SUBROUTINE  SNRcalc(KTsnr^nr,nyquist,n) 

c  THIS  S/R  CALCULATES  THE  AVERAGE  SNR  AS  A  FUNCTION  OF 
c  RADIUS 

REAL*8  KTsiir(n,n),  rsnr(n/2) 

INTEGER! 
n2  =  ii/2 
n2pl  =  ii2  +  1 
DO  10  r  =  1,  nyquist 
nsnr  =  0 
DO  20  i  =  1,  n 
DO20j  -  l,n 
X  =  fioat(j  -  (n2pl)) 
y  =  fioat(i  -  (n2pl)) 
radius  =  sqrt(x**2.0  +  y**2.0) 

IF  ((radius.GT.float(r-l))AND. 

+  (radius.LE.float(r)))  THEN 
nsnr  =  nsnr  +  1 
rsnr(r)  =  rsnr(r)  +  KTsnr(iJ) 

ENDIF 

20  CONTINUE 

rsnr(r)  =  rsnr(r)/nsnr 
10  CONTINUE 
RETURN 
END 
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SUBROUTINE  Star(DATA,n) 

c  THIS  S/R  CREATES  A  SIMULATED  POINT  SOURCE  OR  STAR 

COMPLEX*16  DATA(n,n) 
n2pl  =  i)/2  +  1 
DO  10  i  =  1,  n 
DO  10  j  =  l,n 
X  =  floatQ  -  n2pl) 
y  =  float(n2pl  -  i) 

IF  ((x.EQ.0.0)  AND.(y.EQ.0.0))  THEN 
DATA(id)  =  (1.0, 0.0) 

ELSE 

DATA(ij)  =  (0.0, 0.0) 

ENDIF 

10  CONTINUE 
RETURN 
END 
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SUBROUTINE  Writeffle(DATAl,data2,fflel,file2, 

+  nyquistyn) 

c  THIS  S/R  WRITES  THE  DATA  TO  A  FILE  IN  THE  FORMAT  REQUIRED 

c  BY  A  PROGRAM  WHICH  PRODUCES  A  CONTOUR  PLOT  OF  THE 
c  IMAGE 

COMPLEX*16  DATAl(n,n) 

REAL*8  data2(n/2-l) 

REAL  X,  lambda,  RO,  pi 
INTEGER  r 

CHARACIER*16  filel,  ffle2 

COMMON  A^ARS2/DIAM,OBSCUR,LAMBDA,RO,SECDIM 
PIXSCALE,TPFDIM, FILENAME 

OPEN  (UNIT=30,FILE=filel,STATUS=’NEW’) 

OPEN  (UNIT=40,FILE=file2,STATUS=’NEW’) 
pi  =  acos(-1.0E+00) 

DO  10  i  =  1,  n 
DO  10  j  =  l,n 
WRnE(30,*)  DATAl(ij) 

10  CONTINUE 

DO  20  r  *  1,  nyquist 

X  =  r  *  (lamb^pixscale)  *  (180/pi)  *  3600.0 
WRrrE(40,*)  X,  data2(r) 

20  CONTINUE 
RETURN 
END 
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APPENDK  E.  SPECIFIC  KNOX-THOMPSON  SUBROUTINES 


c  THE  FOLLOWING  SUBROUTINES  WERE  GENERATED  BY  THE 
c  THESIS  AUTHOR  AND  ARE  USED  BY  THE  KNOX-THOMPSON 
c  PROGRAM  ONLY. 


SUBROUTINE  UST 


SUBROUTINE  Initialize(OBJDATA,GAUSSIAN,PUPIL,F, 

+  TEMPDATA,filel,file2,fwd4nv, 

+  lii,iiframes,nhoton,i^quist, 

+  ofEset,cent,n) 

c  THIS  S/R  INITIALIZES  SOME  OF  THE  PROGRAM  VARIABLES  BY 
c  QUEARYING  THE  USER  FOR  INPUT.  IT  ALSO  CREATES  THE 
c  OBJECT  DESIRED  BY  THE  USER  BY  CALLING  EITHER  ASTEROID, 
c  BISTAR,  OR  STAR  S/R. 

INTEGER  ofEset 
CHARACrER*16  filel,  ffle2 
CHARACTER*!  cent,  object 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  INTEGER  OFFSET  (<=  4  AND  <=  #’ 
WRITE(*,*)TIXELS/RO):’ 

READ(*,*)of6set 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  INTEGER  NUMBER  OF  SNAPSHOTS:’ 

READ(*,*)n£rames 

WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  INTEGER  NUMBER  OF  OBJECT  PHOTONS  PER’ 
WRrTE(*,*)’SNAPSHOT:’ 

READ(*,*)nphoton 
WRTTEC*,*)’  ’ 

WRrTE(*,*)’ENTER  INTEGER  NYQUIST  VALUE:’ 

READ(*,*)nyquist 
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10  WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  OBJECT  TYPE:  "  A  "  FOR  ASTEROID,’ 
Write(*,*)"'  b  "  FOR  BINARY  STAR,  OR  "  S  "  FOR  STAR.’ 
READ(*,40)  object 

IF  ((object.EQ.’a’).OR.(object.EQ.’A’))  THEN 
CALL  AST(OBJDATA,GAUSSIAN,PUPIL,F,TEMPDATA. 

+  fwd,iiiv4n,ii) 

ELSEIF  ((objcct.EQ.’b’).OR.(object.EQ.’B’))  THEN 
CALL  Bistar(OBJDATA,n) 

ELSEIF  ((object.EQ.’s’).OR.(object.EQ.’S’))  THEN 
CALL  Star(OBJDATA,n) 

ELSE 

WRITE(*,*)’INCORRECT,  REENTER’ 

GOTO  10 
ENDIF 

20  WRITE(*,*)’  ’ 

^'RITE(*,*)’WOULD  YOU  LIKE  THE  DATA  CENTROIDED?  (Y/N)’ 
READ(*,40)  cent 

IF  ((cent.NE.’y’)j\ND.(centNE.’Y’)AND. 

+  (cent.NE.’n’).AND.(cent.NE.’N’))  THEN 
WRITE(*,*)’INCORRECT,  REENTER’ 

GOTO  20 
ENDIF 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  KT  OUTPUT  FILE  NAME  (16  CHAR’ 
WRITE(*,*)’MAX):’ 

READ(*,30)  filel 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  KT  SNR  FILE  NAME  (16  CHAR’ 
WRriE(*,*)’MAX):’ 

READ(*,30)  ffle2 
30  FORMAT(A16) 

40  FORMAT(Al) 

RETURN 

END 
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+  + 


SUBROUTINE  KTrecon(KTPHASOR,I,IKT,KTsnr, offset, 
xvarTkt,xvarikt,icounter, 
nframes,nyquist,n) 

c  THIS  S/R  IS  THE  HEART  OF  THE  KT  RECONSTRUCTION  PROGRAM, 

c  IT  DETERMINES  THE  VARIANCE-WEIGHTED,  NOISE-BIAS- 
c  CORRECHED  (TIOSS-SPECTRUM  OF  THE  DEGRADED  IMAGE  AND 
c  RECONSTRUCTS  THE  PHASEORS  FROM  THE  AVERAGE  OF  THIS 
c  CROSS-SPECTRUM.  THIS  S/R  IS  CALLED  ONCE  FOR  EACH 
c  OBJECT/POINT  SOURCE  SNAP  SHOT  TO  DETERMINE  A  RUNNING 
c  AVERAGE  CTIOSS-SPECTRUM  AND  WITH  THE  LAST  SNAP  SHOT, 
c  THE  PHASOR  IS  RECONSTRUCTED. 

COMPLEX*16  KTPHASOR(-n/2:ii/2-l,-n/2:n/2-l),  CHEMP, 

+  IKT(-n/2:n/2-l,-n/2:n/2-l,-0:4,-4:4), 

+  I(-n/2:n/2-l,-n/2:n/2-l),  PTEMPl,  PTEMP2 
REAL*8  KTsnr(-n/2:ii/2-l,-n/2:n/2-l),  sigmar,  sigmai, 

+  xvarrlct(-n/2:ii/2-l,-ii/2:n/2-l,-0:4,-4:4),  snr, 

+  xvarikt(-ii/2:ii/2-l,-n/2:n/2-l,-0:4,-4:4),  sigma, 

INTEGER  r,  di,  dj,  offset 
k  =  0 

KTPHASOR(k,k)  =  (1.0,0.0) 

DO  10  r  =  1,  nyquist 
DO  10  u  =  0,  r 
IF  (u.EQ.0)  THEN 
lim  =  0 
ELSE 
lim  =  -r 
ENDIF 

DO  10  jj  =  lim,  r 

rad  =  sqrt(float(ii)**2.0  +  float(jj)**2.0) 

IF  ((rad.LE.float(r)).AND.(rad.GT.float(r-l)))  THEN 
IF  (icounter.EQ.n£rames)  nsnr  =  0 
DO  20  di  =  0,  offset 
DO  20  dj  =  -offset,  offset 
drad  =  sqrt(float(di)**Z0  +  float(dj)**2.0) 

IF  ((drad.LE.float(r))AND. 

-I-  (drad.LE.float(offset)))  THEN 
idi  =  ii  -  di 
jdj  =  jj  -  dj 

radd  =  sqrt(float(idi)**2.0  +  float(jdj)**2.0) 
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++  ++  ++  ++  ++ 


IF  (raddXE.float(r-l))  THEN 
iHC  =  -ii 
jHC  =  -jj 

IF  (I(idijdj).EQ.(0.0,0.0)) 

+  I(idijdj)  =  (1.0, 0.0) 

IF  (I(iijj).EQ.(0.0,0.0)) 

+  I(iijj)  =  (1.0, 0.0) 

CTEMP  =  (I(idiodj)  •  DCONJG(I(iyi)))  - 
+  DCX)NJG(I(di,dj)) 

IKT(idijdj,di,dj)  =  IKT(idydj,di,dj)  + 

+  CTEMP 

xvarrkt(idijdj,di,dj)  = 
xvarTkt(idijdj,di,dj)  + 

(DREAL(CrEMP))**10 
xvarikt(idijdj,di,dj)  = 
xvarikt(idydj,di,dj)  + 

(DIMAG(CTEMP))**2.0 
IF  (icounter.EQ.nframes)  THEN 
nsnr  =  nsnr  +  1 

sigmar  =  dabs((xvarrkt(idijdj,di,dj)  - 
((DREAL(IKT(idijdj,di,dj)))**2.0)/ 
(^oat(n£rames))/dfloat(n£t^es  •  1)) 
sigmai  -  dabs((xvarikt(idijdj,di,dj)  • 
((DIMAG(IKT(idiJdj,di,dj)))**2.0)/ 
dfloat(iijErames))/dfloat(nframes  •  1)) 
sigma  =  dsqrt(sigmar  +  sigmai) 
snr  =  (((ABS(IKT(idiodj,di,dj)))/ 
dfloat(n£rames))/sigma)  * 
dsqrt(dfloat(ii&^es)) 

KTsnr(iijj)  =  KTsnr(iijj)  +  snr 
PTEMPl  =  IKT(idijdj,di,dj)/ 

+  ABS(IKT(idijdj,di,dj)) 

PTEMP2  =  KTPHASOR(idijdj) 
KTPHASOR(ujj)  =  KTPHASOR(iijj)  + 

+  ((snr**2.0)  *  (DCONJG(PTEMPl/PTEMP2))) 

ENDIF 
ENDIF 
ENDIF 

20  CONTINUE 
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IF  (icounter.EQ.n£rames)  THEN 
ICTsnr(iiji)  =  KTsnr(iijj)/iisnr 
KTsnr(iHCjHC)  =  KTsnr(iijj)/nsnr 
KTPHASOR(iiaj)  =  KTPHASOR(iiaj)/ 

+  ABS(KTPHASOR(iijj)) 

ICrPHASOR(iHCjHC)  =  DCONJG(KTPHASOR(iijj)) 
ENDIF 
ENDDF 

10  CONTINUE 
RETURN 
END 
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APPENDK  F.  SPECmC  TRIPLE^ORREIATION  SUBROUTINES 


c  THE  FOLLOWING  SUBROUTINES  WERE  GENERATED  BY  THE 
c  THESIS  AUTHOR  AND  ARE  USED  BY  THE  TRIPLE-CORRELATION 
c  PROGRAM  ONLY. 


SUBROUTINE  LIST 


SUBROUTINE  BSrecon(BSPHASOR,I,IBS,IDBS,BSsnr,o£&et, 

-I-  xvarrbs,xvaribs,icounter,nframes, 

+  i^quist,nphoton,n) 

c  THIS  S/R  IS  THE  HEART  OF  THE  BS  RECONSTRUCTION  PROGRAM, 

c  rr  DETERMINES  THE  VARIANCE-WEIGHTED, 
c  NOISE-BIAS-CORRECTED BISPECTRUM  OFTHE DEGRADED  IMAGE 

c  AND  RECONSTRUCTS  THE  PHASEORS  FROM  THE  AVERAGE  OF 
c  THIS  BISPECTRUM  THIS  S/R  IS  CALLED  ONCE  FOR  EACH 
c  OBJECT/POINT  SOURCE  SNAP  SHOT  TO  DETERMINE  A  RUNNING 
c  AVERAGE  BISPECTRUM  AND  WITH  THE  LAST  SNAP  SHOT,  THE 
c  PHASOR  IS  RECONSTRUCTED. 

COMPLEX*16  BSPHASOR(-n/2:n/2-l,-ii/2;n/2-l),  CTEMP, 

+  PTEMP2,  IBS(-n/2:n/2-l,-n/2:n/2-l,0:4,-4:4), 

+  PTEMPl,  I(-n/2:n/2-l,-n/2:n/2-l),  ID^0:4,-4:4) 

REAL*8  BSsnr(-n/2:n/2-l,-n/2:n/2-l),  sigmar,  sigmai, 

+  snr,  xvarrbs(-ii/2:n/2-l,-n/2:n/2-l,0:4,-4:4), 

+  sigma,  xvaribs(-n/2:ii/2-l,-n/2:ii/2-l,0:4,-4:4) 

INTEGER  r,  di,  dj,  offset 
k  =  0 
loop  =  1 

BSPHASOR(k,k)  =  (1.0,0.0) 

DO  10  r  =  1,  nyquist 
DO  10  ii  =  0,  r 
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+  +  +  + 


IF  (ii.EQ.O)  THEN 
lim  s=  0 
ELSE 
lim  =  -r 
ENDIF 

DO  10  jj  =  lim,  r 

rad  =  sqrt(float(ii)**2.0  +  float(ij)**2.0) 

IF  ((rad.LE.float(r))j\ND.(rad.GT.float(r-l)))  THEN 
IF  (icounter.EQ.iifiames)  nsnr  =  0 
DO  20  di  =  0,  of6set 
DO  20  dj  =  -offset,  offset 
IF  (loop.EQ.l)  THEN 
IDBS(di,dj)  =  IDBS(di,dj)  +  I(di,dj) 

ENDIF 

drad  =  sqrt(float(di)**2.0  +  float(dj)**2.0) 

IF  ((drad.LE.float(r))AND. 

+  (drad.LE.float(offset)))  THEN 
idi  =  ii  -  di 
jdj=jj-dj 

radd  =  sqrt(float(idi)**2.0  +  float(jdj)**2.0) 

IF  (radd.LE.float(r-l))  THEN 
iHC  =  -ii 
jHC  =  -jj 

IF  (I(ididdj).EQ.(0.0,0.0)) 

+  I(ididdj)  =  (1.0, 0.0) 

IF  (I(iijj).EQ.(0.0,0.0)) 

+  I(iijj)  =  (1.0,0.0) 

CTEMP  =  I(idijdj)  * 

+  DCONJG(I(iijj))  • 

+  I(di,dj)  -  (ABS(I(ididdj)))**2.0  - 

+  (ABS(I(iijj)))**2.0  - 

+  (ABS(I(di,dj)))**2.0  +  (2.0  *  nphoton) 

IBS(idiJdj,di,dj)  =  IBS(ididdj,di,dj)  -f- 
+  CTEMP 

xvaiTbs(idijdj,di,dj)  = 
xvarrbs(ididdj,di,dj)  + 
(DREAL(CTEMP))**2.0 
xvaribs(idijdj,di,dj)  = 
xvaribs(idijdj,di,dj)  + 
(DIMAG(CrEMP))**2.0 
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+  +  +  +  +  + 


IF  (icounter.EQ.nframes)  THEN 
nsnr  =  nsnr  +  1 

sigmar  =  dabs((xvarrbs(idijdj,di,dj)  - 
((DREA4IBS(idiodj,di,dj)))**2.0)/ 
dfloat(nframes))/dfloat(n£r^es  -  1)) 
sigmai  -  dabs((xvaribs(idydj,di,dj)  - 
((DIMAG(IBS(idijdj,di,dj)))**2.0)/ 
dfloat(nfiraines))/dfloat(n£ranies  - 1)) 
sigma  -  dsqrt(sigmar  +  sigmai) 
snr  =  (((ABS(IBS(idijdj,di,dj)))/ 
d£loat(nfimnes))/sigma)  * 
dsqrt(dfloat(ii£rames)) 

BSsnr(iijj)  =  BSsnr(iyj)  +  snr 
PTEMPl  =  IBS(idiJdj,di,dj)/ 

+  ABS(IBS(idgdj,di,dj)) 

PTEMP2  =  BSPHASOR(idi  jdj) 

+  *  IDBS(di,dj)/ABS(IDB^di,dj)) 

BSPHASOR(ujj)  =  BSPHASOR(iijj)  + 

+  ((siir**2.0)  *  (DCONJG(PTEMPl/PTEMP2))) 

ENDIF 
ENDIF 
ENDIF 

20  CONTINUE 

IF  (loop.EQ.l)  loop  *  2 
IF  (icounter.EQ.n£rames)  THEN 
BSsnr(iijj)  =  BSsnr(iijj)/nsnr 
BSsnr(iHCjHC)  s  BSsnr(iijj)/nsnr 
BSPHASOR(iijj)  =  BSPHASOR(iyj)/ 

+  ABS(BSPHASOR(uJj)) 

BSPHASOR(iHCdHC)  =  DCONJG(BSPHASOR(iijj)) 
ENDIF 
ENDIF 

10  CONTINUE 
RETURN 
END 
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SUBROUTINE  Initialize(OBJDATA,GAUSSIAN,PUPIL,F, 

+  TEMPDATA,£ilel^e2,fwd,iiiv,lii, 

+  n&’ames,nphoton,nyqiiist,offset, 

+  cent,n) 

c  THIS  S/R  INITIALIZES  SOME  OF  THE  PROGRAM  VARIABLES  BY 
c  QUEARYING  THE  USER  FOR  INPUT.  IT  ALSO  CREATES  THE 
c  OBJECT  DESIRED  BY  THE  USER  BY  CALLING  EITHER  ASTEROID, 
c  BISTAR,  OR  STAR  S/R. 

INTEGER  offset 
CHARACTER*16  filel,  file2 
CHARACTER*!  cent,  object 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  INTEGER  OFFSET  (<=  4):’ 

READ(*,*)offset 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  INTEGER  NUMBER  OF  SNAPSHOTS:’ 

READ(*,*)n£rames 

WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  INTEGER  NUMBER  OF  OBJECT  PHOTONS  PER’ 
WRITE(*,*)’SNAPSHOT:’ 

READ(*,*)nphoton 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  INTEGER  NYQUIST  VALUE:’ 

READ(*,*)nyquist 
10  WRITE(*,*)’  ’ 

WRnE(*,*)’ENTER  OBJECT  TYPE:  "  A  "  FOR  ASTEROID,’ 
WRITEC*,*)"'  B  "  FOR  BINARY  STAR,  OR  "  S  "  FOR  STAR.’ 

READ(*,40)  object 

IF  ((object.EQ.’a’).OR.(object.EQ.’A’))  THEN 

CALL  ast(objdatagaussian,pupil,f,tempdata 

+  fwd,mv,ln,n) 

ELSEIF  ((object.EQ.’b’).OR.(object.EQ.’B’))  THEN 
CALL  Bistar(OBJDATAn) 

ELSEIF  ((object.EQ.’s’).OR.(object.EQ.’S’))  THEN 
CALL  Star(OBJDATA,n) 

ELSE 

WRrTE(*,*)’INCORRECT,  REENTER’ 

GOTO  10 
ENDIF 
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20  WRITEC,*)’  ’ 

WRITE(*,*)’WOULD  YOU  LIKE  THE  DATA  CENTROIDED?  (Y/N)’ 
READ(*,40)  cent 

IF  ((ccnt.NE.y)AND.(ccnt.NE.*Y’)AND. 

+  (ccnt.NE.’n’)j\ND.(ccnt.NE.’N’))  THEN 
WRITE(*,*)’INCORRECT,  REENTER’ 

GOTO  20 
ENDIF 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  TC  OUTPUT  FILE  NAME  (16’ 
WRITE(*,*)’(CHAR  MAX):’ 

READ(*,30)  filel 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  TC  SNR  FILE  NAME  (16  CHAR’ 
WRITE(*,*)’MAX):’ 

READ(*,30)  file2 
30  FORMAT(A16) 

40  FORMAT(Al) 

RETURf4 

END 
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+  + 


APPENDIX  G.  ATMOSPHERIC  DEGRADATION  SUBROUTINES 


c  THE  FOLLOWING  SUBROUTINES  WERE  PROVIDED  BY  CAPTAIN 
c  CHUCK  MATSON  AND  MS.  IDA  DRUNZER  OF  WL/ARCI  AND 

c  MODIFIED  FOR  USE  WITH  THE  KNOX-THOMPSON  AND 
c  TRIPLE-CORRELATION  IMAGE  RECONSTRUCTION  PROGRAMS, 
c  THESE  SUBROUTINES  ARE  REQUIRED  FOR  BOTH  PROGRAMS. 


SUBROUTINE  LIST 


SUBROUTINE  PHSUB(SEIMG,F,TEMPDATA,OBJARR, 
icounter,nframes,nyquist, 
fwd,inv,b,n,mseed) 

c  FUNCTION:  Creates  an  image  of  a  object  located  in  space  by  filtering  a 
c  single  snapshot  of  an  object  with  a  phase  screen. 

c  ORDER  OF  EVENTS  TO  OBTAIN  IMAGE: 

c  -  Create  an  object 

c  -  Create  Gaussian  Array 

c  -  Expose  array  to  what  the  atmosphere  will  do  to  the  object 

c  (Correlation  Filter  Function) 

c  -  Expose  array  to  the  stipulations  of  the  telescope 

c  -  Autocorrelate  the  object  (incoherent  transfer  function) 

c  -  Multiply  the  object  by  the  phase  screen  to  obtain  the  Image 

c  DICTIONARY: 

c  DIAM 
c  ICOUNTER 
c 

c  ISEED 
c 

c  LAMBDA 
c  MU 


-  Telescope  pupil  diameter 

-  Number  of  user  chosen  iterations  of 

phase  screens  to  process 

-  Seed  to  be  used  for  random  number 

generator 

-  Center  wavelength 

-  Mean 
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PDCSCALE  -  Corresponds  to  each  pixel  in  (^cles/m 
RO  -  Variable  characterizing  the  strength  of  atmospheric  turbulence 

with  respect  to  the  optical  system 
TPFDIM  -  Dimension  for  pupil  of  telescope 
PCOUNT  -  The  actual  number  of  photons  per  frame 
NFRAMES  -  The  number  of  frames  to  process 
NPHOTON  -  The  number  of  photons  the  user  wants  per  frame 
NKL  -  The  number  of  Karhunen-Loeve  functions  projected  off  and 

added  back  in  the  phase  map. 

PARAMETER  (IDIM  =  64,LDIM  =  IDIM*IDIM) 

PARAMETER  (IDIM2  =  64,LDIM2  =  IDIM2*IDIM2) 

PARAMETER  (IDIMl  =  .7071*IDIM2) 

PARAMETER  (IWDIM  =  5*IDIM/2,IWDIM2  =  5*IDIM2/2) 

PARAMETER  (N1  =  11, N2  =  (Nl/2)  +  1, 

+  N3  =  (N1  -  1)*(N1  +  2)/2) 

PARAMETER  (N4  =  (N1  -  1)*(N1  +  3)/4) 

PARAMETER  (NKL  =  20) 

CHARACTER*16  FILENAME 
CHARACTER*!  TILT 

COMPLEX*16  PHAMAP(LDIM2),PS1(LDIM),IMAGE(LDIM), 

+  OBJARR(LDIM),PS2(LDIM2),SEIMG(LDIM)  . 

REAL  POIARRAY(LDIM),RVAR,MULTl,MULT2,LAMBDA, 

+  PIXSCALE,ZERO,SIZE,DIAM, TPFDIM, R0,INNERSCALE, 

+  CORR(LDIM2),SECDIM,OUTERSCALE,OBSCUR,RARRAY(LDIM) 

INTEGER  PCOUNT 

DIMENSION  ZKLMAP(IDIM2,IDIM2,NKL),LINPTR(N3),BVAL(N3) 

COMMON  A^ARS/ZERO,RMAX,RMIN,CMINUSl,RCONSTl,FCONSTl,PI 
COMMON  A^ARS2/DIAM,OBSCUR,LAMBDA,RO,SECDIM, 

+  PDCSCALE,TPFDIM,FILENAME 

c  Initialize  variables 

CMINUS1=  (-l.E0,0.E0) 

ZERO  =  0.E0 
MU  =  0.E0 
SIGMA  =  1.E0 


c 

c 

c 

c 

c 

c 

c 

c 
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PI  =  DACOS(-1.0D+00) 

JDIM  =  IDIM 
JDIM2  =  IDIM2 
IC  =  IDIM/2+1 
INNERSCALE  =  .0001 
OUTERSCALE=  10000.E0 
MULTI  =  5.92/INNERSCALE 
MULT2  =  (2.0*PI)/OUTERSCALE 

c  Only  do  certain  subroutines  the  first  time  calling  this  subroutine  (when  wanting 

c  several  snapshots). 

11  CX)NTINUE 

IF((ICOUNTER.EQ.l)AND.(MSEED£Q.l))  THEN 
FILENAME  =  ’phvarld’ 

OPEN(10,FILE=FILENAME,STATUS=’OLD’,ERR=500) 

READ(10,*)DIAM,OBSCUR,LAMBDA,R0 

CLOSE(UNIT=10) 

c  Call  the  subroutine  which  prompts  the  user  for  all  names  and  variables  that  the 

c  subroutine  will  need  in  order  to  run. 

CALL  PARAM2(IDIM,NYQUIST,TILT,ISEED) 

c  Generate  a  filter  function  to  be  used  as  a  multiplier  on  a  complex  gaussian 

c  array  to  be  stored  in  array  CORR.  However,  only  call  the  filter  function  on 

c  multiple  runs  if  there  are  changing  A^ue^  for  the  width  of  the  object  and  RO. 

CALIFILT2(RARRAY,CORR,LDIM2,IDIM2,MULTl,MULT2,PIXSCALE, 
+  ZERO,RO) 

c  Calculate  the  simulated  size  of  the  telescope  (tpfdim).  Compare  it  with  the 
c  dimension  that  is  know  (size).  If  size  is  less  than  tpfdim,  there  will  be  an  error 

c  due  to  the  array  being  too  small  to  store  the  data  from  the  fourier  transform, 

c  and  thus,  data  will  be  truncated. 

SIZE  =  (IDIM/2)-l 
IF(SIZE.LT.TPFDIM)  THEN 
WRITE(*,*)’  ’ 

IF(SIZE.LE,(.5*TPFDIM))  THEN 
WRnE(*,*)’ERROR.  IMAGE  ARRAY  SIZE  MUST  BE 
+  512x512.’ 


143 


ELSE 

WRrrE(V)’ERROR.  IMAGE  ARRAY  SIZE  MUST  BE 
+  256x256.’ 


ENDEF 


WRITE(*,*)’  ’ 


GOTO  400 
ENDIF 


c  Call  BGSCREEN  subroutine  first  time  through  phase  screen.  The  BGSCREEN 
c  subroutine  manipulates  the  phase  map  to  result  in  a  phase  map  with  proper 
c  statistics. 


CALL  BGSCR2(BVAL,ZKLMAP,LINPTR,IDIM2) 
c  End  the  if  statement  if  first  time  throu^  the  subroutine, 
mseed  =  mseed  +  1 


ENDIF 

c  Get  Gaussian  numbers  and  then  multipty  these  two  numbers  by  the  correlation 
c  filter  function  to  obtain  the  phase  map  (PHAMAP). 

DO  20  J  =  1,LDIM2 

CALL  GAUSS(MU,SIGMA,R1,R24>I,1SEED) 

PHAMAP(J)  =  DCMPLX(Rl*CORR(J),R2*CORR(J)) 

20  CONTINUE 

c  Get  phase  screen  which  is  in  the  frequency  domain.  Nrad  perform  the  inverse 

c  FFT  on  array. 

CALL  FFT2D(PHAMAP,TEMPDATA,F,ln,inv,n) 

DO  25  J  =  1,  LDIM 

PHAMAP(J)  =  PHAMAP(J)*(DFLOAT(LDIM)) 

25  CONTINUE 

CALL  PROJ(PHAMAP3VAL,ZKLMAP,LINPTR,iseed, 

+  NKL,R0,N1,N2,N3,N4,IDIM2^I,PDCSCALE,TILT) 
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c  Calculate  ae**(i*realization  of  phase  screen).  Take  cosine  and  sine  of  eveiy 
c  real  element  in  the  array  and  put  into  the  phase  screen  array. 

DO  60  J  =  1  JX)IM2 

PS2(J)  =  DCMPLX(DCOS(DREAL(PHAMAP(J))), 

+  DSIN(DREAL(PHAMAP(J)))) 

60  CONTINUE 

c  Perform  the  Telescope  Pupil  Function  (TPF)  on  the  phase  map  to  obtain  the 
c  coherent  transfer  function.  First,  call  this  routine  to  cut  out  chunk  of  phase 

c  screen  based  on  Telescope  Pupil  Function  of  D/LAMBDA  to  be  stored  in  array 

c  PSl. 

CALL  TPF2A(PS14*S2,IDIM,IDIM2,TPFDIM,SECDIM) 

c  Calculate  the  incoherent  transfer  function  by  multipfying  phase  screen  by  the 
c  auto  correlation  filter  function.  Begin  by  taking  the  Fourier  Transform  of  the 
c  phase  screen. 

CALL  FFT2D(PSl,TEMPDATA,F,ln,fwd,n) 

c  Take  the  magnitude  squared  of  the  array  and  put  the  data  into  the  real  part  of 
c  the  phase  screen. 

DO  120  J  =  1,LDIM 

PSl(J)  =  DCMPLX(DREAL(DCONJG(PSl(J))*PSl(J)),0.e0) 

120  CONTINUE 

CALL  FFT2D(PSl,TEMPDATA,F,ln,inv,n) 

c  Divide  the  phase  screen  by  the  value  at  dc  in  order  to  normalize  the  array  by 
c  that  value  (dcval). 

TEMP  =  DREAL(PS1(IC+(IC-1)*IDIM)) 

DO  125  J  =  1,LDIM 
PS1(J)  =  PS1(J)/TEMP 
125  CONTINUE 
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c  Incoherent  transfer  function  is  completed.  Now  get  the  image  in  the  fourier 

c  domain  multiplying  the  object  and  the  phase  screen  to  get  the  image. 

DO  160  J  =  IJJDIM 
IMAGE(J)  =  (PSl(J)*OBJARR(J)) 

160  CONTINUE 

CALL  FFr2D(IMAGE,TEMPDATA4?,ln4nv^) 

c  Call  the  function,  poisson,  which  will  return  as  a  floating  point  number  in 
c  integer  value  that  is  a  random  deviate  drawn  from  a  Poisson  distribution  of 
c  mean  equal  to  image  value,  using  another  real  function,  uniform,  as  a  source  • 

c  of  uniform  random  deviates. 

PCOUNT  =  ZERO 

DO  230  J  =  1,LDIM 
RVAR  =  DREAL(IMAGE(J)) 

IF(RVAR.LE.0.0)  THEN 
POIARRAY(J)  =  ZERO 

FT 

ITEMP  =  POISSON(RVAR,ISEED) 

POIARRAY(J)  =  DFLOAT(ITEMP) 

PCOUNT  «  PCOUNT  +  ITEMP 
ENDIF 

SEIMG(J)  =  DCMPLX(POIARRAY(J),ZERO) 

230  CONTINUE 

300  RETURN 
400  WRITE(V)” 

WRnE(*,*)’STOPPING  PROGRAM  DUE  TO  ERROR’ 

STOP 

c  Error  Statments 

500  CONTINUE 
WRITEC,*)’  ’ 

WRITE(V)’ERROR,  THIS  FILE  DOES  NOT  EXIST.  REENTER.’ 

WRITE(*,*)’  ’ 

GOTO  11 
800  CONTINUE 
END 
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c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


90 

95 


SUBROUTINE  BGSCR2(bval,2klmap,linptr,nsid) 

This  program  does  the  things  that  result  in  a  phase  screen  with  proper  statistics 


parameter(lu  =  12,lu4  =  18) 
parameter(nr=  1004ikl=20) 

PARAMETER  (N1  =  11,  N2  =  (Nl/2)  +  1, 
+  N3  =  (N1  -  1)*(N1  +  1)11) 

PARAMETER  (N4  =  (N1  -  1)*(N1  +  3)/4) 
real  bval(n3),bvec(n2,n3),rmap(nr,n4),nsidl 
dimension  zidmap(nsid,nsid,nbd) 
dimension  ival(n3,3),ir(n3),linptr(n3) 
character*  16  foame, noname 
nsidl  =  .7071*(nsid  - 1) 
hiame  =  ’eigen.d’ 
noname  =  Idrad-d* 


Dictionary 


bval 

bvec 

ir 

ival 

linptr 

nl 

nkl 

nr 

nsid 

rmap 

zklmap 


contains  eigenvalues  of  covariance  matrix 

contains  eigenvectors  of  covariance  matrix 

pointers 

pointers 

pointers 

number  of  azimuthal  orders  to  be  used  in  the  2^mike  functions 
number  of  Karhunen-l^oeve  functions  projected  off  and  added 
back 

number  of  points  in  rmap  (100  is  more  than  enough) 
dimension  of  side  of  screen 

stored  radial  cut  of  Karhunen-Loeve  functions,  computed  in 
michelin 

stored  Karhunen-Loeve  functions 


open(lu,file  =  &iame,form=’formatted’,status=’old’) 
rewind(lu) 
read(lu,*)  linptr 
rcad(lu,*)  bval 
read(lu,*)  ival 
do  95  i=l,n2 
do  90  j=l,n3 
read(lu,*)  bvec(ij) 
continue 
continue 


147 


10  format(6E4.9) 

open(lu4,file»noname,form=’formatted*,status»’old’) 
rewind  (lu4) 
do  80  i=l,nr 
do  85  j=l,n4 
rcad(lu4,*)  rmap(ij) 

85  continue 
80  continue 
read(lu4,*)  ir 

call  kelp(zklmap,nnap,ir,ival,linptr^4i4, 

+  nsid^,nkl) 

return 

end 


SUBROUTINE  FADD(bmap,zklmap,iq,nsidynkl,sum) 

c  bmap  is  the  phase  screen 

complex*  16  bmap(nsid,nsid)^um 
dimension  zklmap(nsid,nsid,n]d) 
data  pi/3.14159265358979/ 

do  10  i  =  l,nsid 
do  20  j  =  l^id 

bmap(ij)  «  bmap(y)  +  sum*zklmap(y,iq) 

20  continue 
10  continue 

return 

end 


SUBROUTINE  nLT2(RARRAY,C»RR,LDIM2,IDIM2,MULTl, 

+  MULT2,PIXSCALEZEROJRO) 

c  Subroutine  Function:  This  function  generates  an  array  representing  the  effects 
c  of  what  the  atmosphere  will  do  to  an  object 

c  when  it  passes  through  it. 

REAL  RARRAY(LDIM2),CORR(LDIM2), MULTI, MULT2,RCONST, 

+  PIXSCALE,ZERO,RO,FCONST,R,ARRAYW 
INTEGER  JOFFSET 
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c  Calculate  the  radial  average  of  an  input  image. 

c  ARRAYS:  CORR(LDIM2)  -  idim2  x  idim2  array 
c  RARRAY(IDIM2)  >  this  is  a  small  work  array 

c  that  needs  to  be  at  least 

c  IDIM2  in  size 

c  other  variables:  IDIM2  -  x  dimensionl 

c  RCONST,RCONST2  -  real  constants 

c 

c  load  CORR  with  X  (J)  indicies  squared 

FCONST  =  .1517 
DO  10  J=1,IDIM2 
JOFFSET=(J-l)*IDIM2 
R=FLOAT(  J-(IDIM2/2+l)  ) 

RCONST=R*R 
DO  15  I=1,IDIM2 
CORR(I+JOFFSET)  =  RCONST 
15  CONTINIJE 
10  CONTINUE 

c  Add  to  LI,  Y  (I)  indicies  squared  (one  column  of  Y  indicies  stored  in  L2) 

DO  20  I=1,IDIM2 
R=FLOAT(  I.(IDIM2/2+l) ) 

RARRAY(I)  =  R*R 
20  CONTINUE 
DO  30  J=1,IDIM2 
JOFFSET=(J-l)*IDIM2 
DO  40  I  =  1,IDIM2 

CORR(JOFFSET  +  I)  =  RARRAY(I)  +  CORR(JOFFSET  +  I) 

40  CONTINUE 
30  CONTINUE 

c  Move  the  scaling  array  to  another  array  in  order  to  do  calculations  for  the  filter 

c  function 

DO  50  I  =  1,LDIM2 
RARRAY(I)  =  CORR(I) 

50  CONTINUE 
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c  Multiply  the  scaling  array  by  the  correlation  filter  function 


RCX)NST  =  (-11^/12.E0) 

DO  60 1  = 

RARRAY(I)  =  (RARRAY(I)+(MULT2*MlJLT2))**RCONST 
60  CONTINUE 

ARRAYW  =  PIXSCALE*IDIM2 
RCONST  =  (R0/ARRAYW)**(-5.E0/6.E0) 

DO  70 1  =  1,LDIM2 

RARRAY(I)  =  (RARRAY(I)*RCONST)*FCONST 
70  CONTINUE 

c  Set  the  middle  point  in  array  to  zero  to  normalize  the  array  by  that  value 

RARRAY(IDIM2*(IDIM2/2)+(IDIM2/2+l))  =  ZERO 

c  Calculate  second  part  of  correlation  filter  function  and  multiply  it  the  first 
c  part  to  complete  the  filter  array. 

RCONST  =  -2.0*(MULT1*MULT1) 

DO  90  1=  1,LDIM2 

CORR(I)  =  (EXP(CORR(I)/RCONST))*RARRAY(I) 

90  CONTINUE 

900  CONTINUE 
RETURN 
END 


SUBROUTINE  GAUSS(MU,SIGMA,RNUMl,RI»UM2,PI,iseed) 
c  PURPOSE: 

c  GET  GAUSSIAN  DISTRIBUTED  RANDOM  NUMBER  WITH  MEAN  MU 

c  AND  STANDARD  DEVIATION  SIGMA 

REAL  MU,SIGMA,Yl,Y2,RNUMl,RNUM2,PI,TWOPI 

TWOPI  =  2.eO*PI 

Yl=  uniform(iseed) 

Y2=  uniform(iseed) 
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c  Calculate  first  random  number 

RNUMl  =  MU  +  SIGMA*SQRT(-2.0*ALOG(Y1))*COS(TWOPI*Y2) 

c  Calcidate  the  second  random  number 

RNUM2  =  MU  +  SIGMA*SQRT(-ZO*AIjOG(Yl))*SIN(TWOPrY2) 

RETURN 

END 


SUBROUTINE  INPROD(bmap;eklmap,iq,nsid,nkl^p) 

complex*  16  bmap(nsid,nsid) 

complex*  16  zap 

dimension  zklmap(nsid,nsid,nkl) 

zap  =  (0.0, 0.0) 
area  =  0.0 
do  10  i  =  l,nsid 
do  20  j  =  l,nsid 

zap  =  zap  +  bmap(ij)*zklmap(ij,iq) 
area  *=  area  +  (  zkhnap(i  j,iq)  )**2 
20  continue 
10  continue 

zap  =  zap/area 

return 

end 
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SUBROUTINE  KELP(zklmap^ap,ir,ival,lmptr, 

c  zklmap  is  the  map  of  the  k-1  functions,  sigh 

dimension  rmap(nr,n4),ir(n3) 
dimension  zkhnap(nsi^nsid,nkl) 
dimension  ival(n3,3)4inptr(n3) 
data  pi/3.14159265358979/ 
icent  =  nsid/2  +  1 
dx  =  1.0/(nsid/2  - 1) 
dr  =  1.0/(nr  - 1) 
do  15  iq  =  l,nkl 
m  *  ival(  linptr(iq),l  ) 
iod  =  ival(  linptr(iq),2  ) 
ipq  =  ir(iq) 
do  10  i  =  l,nsid 
X  =  (i  -  icent)*dx 

do  20  j  s  l^isid 

zkhnap(ij,iq)  =  0.0 
y  =  (j  -  icent)*dx 
r  =  sqrt(x**2  +  y**2) 
th  =  0.0 
if(  r.le.1.0)  then 
if(  r.gt.0.0)  then 
th  =  at£^(y,x) 
if(th.le.0.0)  th  =  2.0*pi  +  th 
else 

th  =  0.0 
endif 

if(r.lt.0.99999)  then 
kl  =  int(r/dr)  +  1 
k2  =  kl  +  1 
rl  =  (kl  -  l)*dr 
r2  =  (k2  -  l)*dr 
coef  1  =  (r2  -  r)/dr 
coef2  =  (r  -  riydr 

top  =  coefl*rmap(kl,ipq)  +  coef2*rmap(k2,ipq) 
else 

top  =  rmap(nr,ipq) 
endif 

if(iod.eq.l)  zz  =  cos(m*th) 
if(iod.eq.2)  zz  =  sin(m*th) 
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zklmap(y,iq)  =  22*top 
endif 
continue 
continue 
continue 
return 
end 


SUBROUTINE  PARAM2(IDIM^QUIST,TILT,ISEED) 

The  function  of  this  subroutine  is  to  ask  the  user  various  questions  about  the 
subroutines  he  or  she  wishes  to  use  and  what  values  he  or  she  wants  to  assign 
to  the  variables  in  the  program. 

REAL  DIAM,Ij\MBDA,R0,PIXELNUM,OBSCUR4*IXSCALE, 

+  secdihtpfdim 

INTEGER  IDIM 
CHARACTER*  16  FILENAME 
CHARACTER*!  ANSWER,FLAG,TILT 
COMMONA^ARS2/DIAM,OBSCUR,LAMBDA,RO,SECDIM, 

+  PIXSCALE,TPFDIMJTLENAME 

WRITE(*,*)’  ’ 

WRITE(*,*)’THE  DEFAULT  VALUES  FOR  THIS  PROGRAM  ARE  AS 
+  FOLLOWS:  ’ 

WRITE(*,*)’  ’ 

WRrTE(*,*)’  TELESCOPE  DIAMETER:  ’,DIAM 

WRITE(*,*)’  CENTER  WAVELENGTH:  LAMBDA 

WRrTE(*,*)’  RO:  ’,R0 

WRrTE(*,*)’  TELESCOPE  SECONDARY  DIAMETER:’, OBSCUR 

CONTINUE 
WRITEC*,*)’  ’ 

WRrTE(*,*)’USE  THE  DEFAULT  VALUES?  (Y/N)’ 

READ(*,14)  FLAG 

IF(((FLAG.NE.’Y’)J^,(FLAG.NE.’y’))j\ND. 

+  ((FLAG.NE.’N’).AND.(FLAG.NE.’n’)))  THEN 
WRrTE(*,*)’  ’ 

WRITE(*,*)TOUR  ANSWER  HAS  TO  BE  EITHER  Y  OR  N. 

+  REENTER’ 

FLAG  =  ” 

GOTO  20 
ENDIF 
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30  CONTINUE 

IF((FLAG.EQ.’N’).OR.(FLAG.EQ.’n’))THEN 
WRITEC*,*)’  ’ 

WRITE(V)’  A:  ALL  VARIABLES’ 

WRITEC*,*)’  D:  TELESCOPE  DIAMETER’ 

WRITE(*,*)’  W:  CENTER  WAVELENGTH  (LAMBDA)’ 
WRnE(*,*)’  R:  RO’ 

WRITE(*,*)’  S:  SECONDARY  DIAMETER’ 

WRITE(*,*)’  ’ 

WRITE(*,*)’WHICH  VALUE  WOULD  YOU  LIKE  TO  CHANGE?’ 
READ(*,14)  ANSWER 

IF(((ANSWER.NE.’A’)AND.(ANSWER.NR’a’))j\ND. 

+  ((ANSWER.NE.’D’)J^ND.(ANSWER.NE.’d’))j\ND. 

+  ((ANSWERJ^.’L’)AND.(ANSWER.NE.’1’))AND. 

+  ((ANSWER.NE.’W’).AND.(ANSWER.NE.’w’))AND. 

+  ((ANSWEILNE.’R’)AND.(ANSWER.NE.’r’))AND. 

+  ((ANSWER.NE.*S’)AND.(ANSWER.NE.’s’)))THEN 
WRITE(*,*)’  ’ 

WRITER, *)’ANSWER  IS  NOT  A  CHOICE.  REENTER’ 

GOTO  30 
ENDIF 
WRnE(*,*)’  ’ 

WRITEC*,*)*!™  CURRENT  VALUES  ARE:  ’ 

WRITE(*  •)’  ’ 

WRITE(*,*)’  TELESCOPE  DIAMETER:  ’,DIAM 

WRITEC,*)’  CENTER  WAVELENGTH:  ’, LAMBDA 

WRnE(*,*)’  RO:  ’,R0 

WRITE(*,*)’  SECONDARY  DIAMETER  ’,OBSCUR 

WRnE(*,*)’  ’ 

WRITE(*,*)’  ’ 

IF((ANSWER.EQ.’A’).OR.(ANSWER.EQ.’a’))THEN 
WRnE(*,*)’INPUT  THE  FOLLOWING’ 

WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  NEW  TELESCOPE  DIAMETER:’ 
READ(*,100)  DIAM 
WRnE(*,*)’  ’ 

WRITE(*,*)’ENTER  NEW  CENTER  WAVELENGTH:’ 
READ(*,100)  LAMBDA 
WRITEC,*)’  ’ 

WRITE(*,*)’ENTER  NEW  RO:’ 

READ(*,100)  RO 
WRITEC*,*)’  ’ 

WRITER, •)’ENTER  NEW  TELESCOPE  SECONDARY  DIAM:’ 


154 


50 


READ(*,100)  OBSCUR 

ELSEIF((ANSWER.EQ.’D’).OR.(ANSWER.EQ.’d’))THEN 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  NEW  TELESCOPE  DIAMETER:’ 
READ(*,100)DIAM 

ELSEIF((ANSWER.EQ.’W’).OR.(ANSWER.EQ.’w’))THEN 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  NEW  CENTER  WAVELENGTH:’ 
READ(*,100)  LAMBDA 

ELSEIF((ANSWER.EQ.’R’).OR.(ANSWER.EQ.’r’))THEN 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  THE  NEW  RO:’ 

READ(*,100)R0 

ELSEIF((ANSWER.EQ.’S’).OR.(ANSWER.EQ.’s’))THEN 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  THE  NEW  SECONDARY  DIAMETER:’ 
READ(*,100)  OBSCUR 
ENDEF 
CONTINUE 


WRrrE(*,*)’  ’ 

WRITER, •)’WOULD  YOU  LIKE  TO  CHANGE  ANY  OTHER  VALUE? 
+  (Y/N)’ 

READ(*,14)ANSWER 

IF(((ANSWER.NE.’Y’).AND.(ANSWER.NE.’y’)).AND. 

+  ((ANSWER.NE.’N’).AND.(ANSWER.NE.’n’)))THEN 


WRITE(*,*)’  ’ 

WRnE(*,*)’YOUR  ANSWER  HAS  TO  BE  EITHER  Y  OR  N. 
+  REENTER’ 

GOTO  50 
ENDIF 

IF((ANSWER.EQ.’Y’).OR.(ANSWER.EQ.’y’))THEN 
GOTO  30 
ENDIF 

OPEN(UNIT=  10,FILE=FILENAME,FORM=’FORMATTED’, 

+  STATUS=’UNKNOWN’) 

REWIND(IO) 

WRITE(10,*)DIAM, OBSCUR, LAMBDA,R0 
CLOSE(UNIT=10) 

ENDIF 


WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  THE  NUMBER  OF  PIXELS  PER  RO:’ 
READ(*,*)PDCELNUM 
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WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  SEED  VALUE:’ 

READ(*,*)ISEED 
PDCSCALE  =  RO/PDCELNUM 
TPFDIM  »  nmt(DIAM*PDCELNUM/RO) 

SECDIM  =  iimt(OBSCUR*PDCELNUM/RO) 

60  CONTINUE 
65  CONTINUE 
WRITE(*,*)’  ’ 

WRITE(*,*)’IX)  YOU  WISH  TO  INCLUDE  X/Y  TILT?  (Y/N)’ 
READ(*,14)TILT 

IF(((TiTJ^.’Y’)j\ND.(TILT.NE.’y’))AND. 

+  ((TILT.NE.’N’)j\ND.(TILT.NE.’n’)))THEN 
WRITE(*,*)’  ’ 

WRITE(*,*)’YOUR  ANSWER  MUST  BE  EITHER  Y  OR  N. 
+  REENTER’ 

GOTO  65 
ENDIF 

c  Check  the  value  for  nyquist  to  assure  it  is  equal  to  tpfdim 
IF(NYQUIST.NE.TPFDIM)  GOTO  999 
RETURN 

c  Format  Statements 

14  FORMAT(Al) 

100  FORMAT(E11.4) 
no  FORMAT(I6) 

c  Error  statements 

999  WRITE(V)” 

WRITE(*,*)’ERROR.  VALUE  FOR  NYQUIST  MUST  EQUAL 
+  TO:’, TPFDIM 

+  ,’  NOT  ’, nyquist,  ’.’ 

WRITE(*,*)’  ’ 

WRirE(*,*)TROGRAM  STOPPING.  CHANGE  NYQUIST.’ 

STOP 

END 
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SUBROUTINE  PROJ(bmap,bval,zklniap,linptr,iseed, 

+  nkl,r0,nl,n2,n3,n4,idiin2,pi,pixscale,tilt) 

c  The  function  of  this  subroutine  is  to  project  off  the  low-order  Karhunen-Loeve 
c  functions  and  then  adds  them  back  with  the  proper  strength. 

c  nl  is  the  number  of  azimuthal  orders,  bmap  is  the  phase  screen 

complex*  16  bmap(idim2,idim2),zap,sum 
dimension  bval(n3),linptr(n3) 
dimension  zklmap(idim2,idii]^nkl) 
real  x,y,yi,pi,twopi,pixscale,drO,rO 
integer  izem,idi^ 
character  tilt 


izem  =  5 

twopi  =  2.e0*pi 

drO  =  (pixscale*idim2)/r0 

DO  10  i  =  l,izem 

c  A  random  size  is  generated  for  each  karhunen  function 

X  =  sqrt(bval(linptr(i)  )  ) 
y  =  O.eO 
yi  =  O.eO 
mu  =  0.e0 
sigma  =  X 

call  gauss(mu,sigma,y,yi,pi,iseed) 

sum  =  cmplx(Y,YI) 

sum  =  sum*(dr0*  *.83333333) 

call  inprod(bmap,zklmap,i,idim2,nkl,zap) 

c  The  function  is  projected  out 

if(tilt.eq.’N’.or.tilt.eq.’n’)  then 
if(i.le.2)sum  =  0.0 
endif 

sum  =  sum  -  zap 

call  fadd(bmap,zklmap,i,idim2,nkl,sum) 

10  CONTINUE 
RETURN 
END 
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SUBROUTINE  TPF2A(PS1,PS2,IDIM,IDIM2,TPFDIM,SECDIM) 

c  Function:  The  function  of  this  subroutine  is  to  cut  out  a  portion  of  array  with 
c  the  dimension  of  TPFDIM  and  put  the  data  into  an  array  of  zeros. 

COMPLEX*16  PS1(IDIM,IDIM),PS2(IDIM2,IDIM2) 

REAL  TPFDIM,SECDIM 
INTEGER  IC,IC2 

IC  =  IDIM/2  +  1 
IC2  =  IDIM2/2  +  1 

c  Zero  out  psl  to  ensure  that  no  reminents  of  previous  phase  screen  is  left  over. 

DO  10  J  =  1,IDIM 
DO  10  K  =  1,IDIM 
PS1(J,K)  =  (0.0, 0.0) 

10  CONTINUE 

c  Divide  the  telescope  pupil  dimension  2  in  order  to  obtain  the  radius  instead 
c  of  the  diameter  of  the  pupil. 

TPFDIM2  «  TPFDIM/2 
SECDIM2  *  SECDIM/2 

IF(SECDIM.EQ.0.0)PS1(IC,IC)  =  PS2(IC2,IC2) 

c  Now  insert  the  phase  screen  into  the  larger  array.  The  smaller  array  is  the  data 

c  that  this  telescope  is  able  to  see  due  to  the  inhibition  of  its  size. 

DO  20  I  =  0,  TPFDIM2 
DO  20  J  *  0,  TPFDIM2 
IF((r*2  +  J**2).LE.TPFDIM2**2)  THEN 
IF((I**2  +  J**2).GT.SECDIM2**2)  THEN 
PS1(IC+I,IC+J)  =  PS2(IC2+I,IC2+J) 

PS1(IC+I,IC-J)  =  PS2(IC2+I,IC2-J) 

PS1(1C-I,IC+J)  =  PS2(1C2-I,IC2+J) 

PS1(IC-I,IC-J)  *  PS2(IC2-I,IC2-J) 

ENDIF 

ENDIF 

20  CONTINUE 
RETURN 
END 
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FUNCTION  GAMMLN(xx) 

real  cof(6),stp,half^one4pf,x,tmp^r 

data  coUtp^6.18009173d0,-86.50532033d0, 

+  24.01409822dO,-1.231739516dO,.120858003d-2, 
+  -J36382d-5,2J0662827465d0/ 
data  half, one, ^^.5d0,l>0d0,5.5d0/ 


X  =  XX  -  one 

tmp  =  X  + 

tmp  =  (x  +  half)*log(tmp)  -  tmp 
ser  =  one 

do  lOj  =  1,6 
X  =  X  +  one 
ser  =  ser  +  cof(j)/x 
10  continue 

gammln  =  tmp  +  log(stp*ser) 

return 

end 


REAL  FUNCTION  POISSON(pmean,  iseed) 
c  mean  of  poisson  distribution 
real  pmean 

c  Seed  for  random  number  generator 
integer  iseed 

c  Summary  of  purpose 

c  Generate  a  random  number  with  a  poisson  distributionof  mean  pmean  (poisson 
c  deviate) 

c  Author 

c  numerical  recipes,  p207,  routine  poidev.for 
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c  Modifications 
c  13  nov,  1987:  pat  fitch 

c  1988:  erik  m  Johansson  -  fixed  integer  conversion  problem  with  maxint 

c  Routines  called 

c  uniform 
c  gammln 

c  (c)  Copyright  1987  the  Regents  of  the  University  of  California.  All  rights 
c  reserv^. 

c  This  software  is  a  result  of  work  performed  at  Lawrence  Livermore  National 
c  Laboratory.  The  United  States  Government  retains  certain  rights  therein. 

real  pi 

parameter  (pi=3.141592654) 

c  maxint  is  the  largest  real  number  which  can  be  converted  to  an  integer  without 
c  resulting  in  an  arithmetic  error  (32  bits) 

real  maxint 

parameter  (maxint  =  .214748352el0) 

real  pexpmean,  oldmean,  t,  em,  sq,  alxm»  y,  gammln, 

+  uiiiform 
integer  times 
data  oldmean  Al./ 
if(pmean.lt.  12.0)then 

c  Use  direct  computation  method 

if(pmean.ne.oldmean)  then 

c  New  mean,  calculate  the  exponential 

oldmean  =  pmean 
pexpmean=  exp(-pmean) 
endif 
em  =  -1.0 
t  =  1.0 

2  em  =  em  +  1.0 
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t  =  t  *  iiniform(iseed) 
if(tg:t.pexpinean)  go  to  2 
else 

if(pmeanjie.oldmean)  then 
oldmean  =  pmean 
sq  =  sqrt(2.0*pmean) 
ato=  alog(pmean) 

c  Natural  log  of  gamma  function  =  gammln 

pexpmean=  pmean*alxm  -  gammln(pmean+l.) 
endif 


times  =  0 

1  y  =  tan(  pi*  uniform(iseed)) 
times  =  times  +  1 
if  (times  .ge.  1000)  then 

write(*,*)’ERROR:  STUCK  IN  IX>OP  IN  POISSON’ 
stop 
endif 

em  =  sq*y+pmean 

if((em  .It.  0.)  .or.  (em  .gt.  maxint))  go  to  1 
em  =  int(em) 

t  «  0.9*(l.+y**2)*exp(em*alxm-gammln(em+l.)- 
+  pexpmean) 
if(uniform(iseed).gtt)  go  to  1 
endif 

poisson  =  em 

return 

end 
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REAL  FUNCTION  UNIFORM(seed) 
integer*4  seed 
c  Summary  of  purpose 

c  implements  a  multiplicative  linear  congruential  generator  generates  random 
c  numbers  uniformly  distributed  on  the  open  interval  0.0  to  1.0  (0  and  1  are  NOT 
c  included) 

c  Author 

c  from  "Random  Number  Generators:  Good  Ones  Are  Hard  to  Find,"  Stephen 
c  K.  Park  and  Keith  W.  Miller,  Communications  of  the  ACM,  October  1988,  Vol 
c  31,  No  10,  pp  1192  -  1201  routine  is  listed  on  p  1195 

c  modifications 

c  4/11/89  erik  m  Johansson  -  modified  code  to  use  less  computational  steps, 
c  Equations  used  are  from  the  article  "Efficient  and  Portable  Combined  Random 

c  Number  Generators,"  Pierre  L’Ecuyer,  Communications  of  the  ACM,  June 
c  1988,  Vol  31,  No  6,  at  the  top  of  p745. 

c  (c)  Copyright  1989  the  Regents  of  the  University  of  California.  All  rights 
c  reserved. 

c  This  software  is  a  resxilt  of  work  performed  at  Lawrence  Livermore  National 
c  Laboratory.  The  United  States  Government  retains  certain  rights  therein. 

integer*4  a,  m,  q,  r 

parameter  (a  =  16807) 
parameter  (m  =  2147483647) 
parameter  (q  =  127773) 
parameter  (r  =  2836) 

real  h 

parameter  (h  =  1.  /  2147483647.) 
integer*4  k 
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c  The  function 


k  =  seed  /  q 

seed  =  a  *  (seed  -  k*q)-k*r 

if  (seed  .It.  0)  seed  =  seed  +  m 

uniform  =  seed  *  h 

return 

end 
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APPENDK  H.  PHASE  ERROR  AND  LOW  PASS  FILTER  PROGRAM 


c  THIS  PROGRAM  WAS  DEVELOPED  BY  THE  THESIS  AUTHOR  AND 
c  DEVELOPS  A  LOW  PASS  FILTER  IN  THE  FREQUENCY 
c  DOMAIN  TO  FILTER  THE  RECONSTRUCTED  OBJECT  DATA  FILE 
c  FOR  BOTH  THE  KNOX-THOMPSON  AND  TRIPLE-CORRELATION 
c  METHODS.  ADDITIONALLY,  THIS  PROGRAM  DETERMINES  THE 
c  AZIMUTHAL  RMS  PHASE  ERROR  FOR  THE  INPUT  FOURIER 
c  SPECTRA. 

c  THE  FOLLOWING  SUBROUTINES  ARE  REQUIRED  FROM 
c  UNIVERSAL  SUBROUTINES  IN  APPENDDC  D: 

c  Complexconj 

c  FFT 

c  FFT2D 

c  NormFFT 

c  Quadswap 


c  AUTHOR: 
c  COMPL.  DATE: 
c  REASON: 
c 

c  GOAL: 
c 
c 
c 


LT  JAMES  M.  LACKEMACHER 
26  OCTOBER  1990 

COMPLETE  REQUIREMENTS  FOR  A  MASTERS 
DEGREE  IN  PHYSICS 
SIMULATE  OBJECT,  DEGRADE  OBJECT, 
RECONSTRUCT  OBJECT  USING  KNOX-THOMPSON 
AND  TRIPLE-CORRELATION  METHODS  FILTER 
THE  DATA  THE  COMPARE  THE  TWO  METHODS. 


PROGRAM  FREQFILTER 

MAIN  PROGRAM  COMPLEX  VARIABLE  LIST 

c  BSDATA  n  X  n  DIM  ARRAY  REPRESENTING  THE  INPUT 

c  TRIPLE-CORRELATION  RECONSTRUCTED  OBJECT  IN 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


BSLP 

F 

KTDATA 

KTLP 

TEMPDATA 

TRUDATA 


FREQUENCY  SPACE 

n  X  n  DIM  ARRAY  REPRESENTING  THE  LOW  PASS 
FILTERED  TRIPLE-CORRELATION  RECONSTRUCTED 
OBJECT  IN  FREQUENCY  SPACE 
n  DIM  ARRAY  USED  IN  THE  FOURIER  TRANSFORM 
n  X  n  DIM  ARRAY  REPRESENTING  THE  INPUT 
KNOX-THOMPSON  RECONSTRUCTED  OBJECT  IN 
FREQUENCY  SPACE 

n  X  n  DIM  ARRAY  REPRESENTING  THE  LOW  PASS 
FILTERED  KNOX-THOMPSON  RECONSTRUCTED 
OBJECT  IN  FREQUENCY  SPACE 
n  X  n  DIM  ARRAY  THAT  IS  USED  AS  A  TEMPORARY 
ARRAY  IN  THE  FOURIER  TRANSFORM 
n  X  n  DIM  ARRAY  REPRESENTING  THE  TRUTH  DATA 


c 


MAIN  PROGRAM  REAL  VARIABLE  LIST 


r 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


bserror 

bslpmod 

bsr 

bssnr 

bstninc 

kterror 

ktlpmod 

ktr 

ktsnr 


n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  SQUARE 
PHASE  ERROR  FOR  THE  TRIPLE-CORRELATION 
RECONSTRUCTED  IMAGE 

n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  MODULUS 
OF  THE  LOW  PASS  FILTERED  TRIPLE-CORRELATION 
RECONSTRUCreD  IMAGE 

n/2  DIM  ARRAY  THAT  REPRESENTS  THE  RADIAL 
FREQ  VALUE  IN  INVERSE  ARCSECONDS  FOR  THE 
TRIPLE-CORRELATION  IMAGE 
n/2  DIM  ARRAY  THAT  REPRESENTS  THE  INPUT  SNR 
VALUES  OF  TRIPLE-CORRELATION  IMAGE 
TRIPLE-CORRELATION  RADIALTRUNCATION  VALUE 
n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  SQUARE 
PHASE  ERROR  FOR  THE  KNOX-THOMPSON 
RECONSTRUCTED  IMAGE 

n  X  n  DIM  ARRAY  THAT  REPRESENTS  THE  MODULUS 
OF  THE  LOW  PASS  FILTERED  KNOX-THOMPSON 
RECONSTRUCTED  IMAGE 

n/2  DIM  ARRAY  THAT  REPRESENTS  THE  RADIAL 
FREQ  VALUE  IN  INVERSE  ARCSECONDS  FOR  THE 
KNOX-THOMPSON  IMAGE 

n/2  DIM  ARRAY  THAT  REPRESENTS  THE  INPUT  SNR 
VALUES  OF  KNOX-THOMPSON  IMAGE 
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c  kttnmc 
c  lambda 
c  rBSerror 
c 
c 

c  rKTeixor 

c 

c 

c  RO 


KNOX-THOMPSON  RADIAL  TRUNCATION  VALUE 
WAVELENGTH  OF  QUASI-MONOCHROMATIC  UGHT 
n/2  DIM  ARRAY  THAT  REPRESENTS  THE  AZIMUTHAL 
RMS  PHASE  ERROR  OF  THE  TRIPLE-CORRELATION 
RECONSTRUCTED  IMAGE 

n/2  DIM  ARRAY  THAT  REPRESENTS  THE  AZIMUTHAL 
RMS  PHASE  ERROR  OF  THE  KNOX-THOMPSON 
RECONSTRUCTED  IMAGE 
COHERENCE  LENGTH 


c  MAIN  PROGRAM  INTEGER  VARIABLE  UST 


c  fwd 
c  iktcount 
c  inv 
c  itccount 
c  In 
c  n 
c  nyquist 
c 
c 

c  numpix 


VALUE  OF  1  FOR  FORWARD  FFT 

NUMBER  OF  KT  PIXELS  WITH  SNR  GREATER  THAN  1.0 

VALUE  OF  -1  FOR  INVERSE  FFT 

NUMBER  OF  TC  PIXELS  WITH  SNR  GREATER  THAN  1.0 

2 '"In  FOR  USE  WITH  FFT  SUBROUTINE 

DIMENSION  OF  ONE  SIDE  OF  2-DIM  ARRAY 

EQUAL  TO  THE  TELESCOPE  PUPIL  FUNCTION 

DERIVED  FROM  THE  FOLLOWING  FORMULA: 

nyquist  =  (telescope  diameter  x  number  of  pixels/r0)/r0 

THE  NUMBER  OF  PIXELS  PER  RO 


c  MAIN  PROGRAM  INTEGER  VARIABLE  LIST 


c 

end 

CHAR  INPUT  TO  DETERMINE  WHETHER  TO  STOP  THE 

c 

PROGRAM 

c 

mean 

CHAR  INPUT  TO  DETERMINE  WHETHER  TO  FIND  THE 

c 

AMSPE 

c 

rite 

CHAR  INPUT  TO  DETERMINE  WHETHER  TO  WRITE 

c 

DATA  TO  FILE 

f 
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MAIN  PROGRAM 


PARAMETER(n=644n=6,£wd=  l,iiiv=-l) 


COMPLEX*  16  KTDATA(n,n),  BSDATA(n,n),  TEMPDATA(n,n), 
+  F(n),  KTLP(n,n),  BSLP(n,n),  TRUDATA(n,n) 

REAL*8  ktlpmod(n,n),  bslpmcxl(n,n),  kterror(n,n), 

+  bseiTor(n,n),  rCTerror(n/2),  rBSeiTor(ii/2), 

+  ktsnr(ii/2),  bssnr(n/2) 

REAL  ktr(n/2),  bsr(i)/2),  kttrunc,  bstninc,  lambda 
CHARACTER*!  mean,  end,  rite 


c  INPUT  THE  DATA 


CALL  Readffle(KTDATA,BSDATA,n) 

WRITE(*,*)’ENTER  NYQUIST  VALUE:’ 

READ(*,*)nyquist 
WRin^*,*)’  ’ 

WRITE(*,*)’ENTER  THE  NUMBER  OF  PDCELS  PER  RO:’ 
READ(*,*)  numpix 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  RO  VALUE  IN  METERS:’ 

READ(*,*)  RO 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  WAVELENGTH  VALUE  IN  METERS:’ 
READ(*,*)  lambda 


WRITE(*,*)’  ’ 

CALL  Readsnr(ktsnr,bssnr,ktr,bsr,nyquist,n) 


c  LOW  PASS  FILTER  BASED  ON  RADIAL  TRUNCATION  VALUE 
c  DETERMINED  FROM  SNR  DATA 


CALL  Tnincval(ktsnr,ktr,kttrunc,nyquist,iktcount,n) 
WRrrE(*,*)’KT  RADIAL  TRUNCATION  VALUE  IS:’ 
WRITE(*,*)  kttnmc,’l/ARCSEC.’ 

WRITE(*,*)’VALUES  GREATER  THAN  THIS  RADIUS’ 
WRITE^,*)’WILL  BE  TRUNCATED.’ 

WRITEC*,*)’  ’ 
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CALL  Truncval(bssnr,bsr,bstrunc,nyquist,itccount,n) 

WRrrE(*,*)TC  RADIAL  TRUNCATION  VALUE  IS:’ 

WRITE(*,*)  bstrunc,’l/ARCSEC’ 

WRITE(*,*)’VALUES  GREATER  THAN  THIS  RADIUS’ 
WRITE(*,*)’WILL  BE  TRUNCATED.’ 

WRITE(* ,*)’  ’ 

c  TRUNCATE  THE  DATA 

CALL  Tnincate(KTDATA,KTLP,kttruiic,lainbda,RO,numpix,n) 
CALLTruncate(BSDATA,BSLP,bstrunc,lambda,RO,numpix,n) 

c  INVERSE  FFT  THE  DATA  AFTER  FILTERING 

CALL  FFT2D(KTLP,TEMPDATA,F,ln,mv,n) 

CALL  FFT2D(BSLP,TEMPDATA,F,ln,inv,n) 

c  DETERMINE  THE  MODULUS  OF  THE  DATA  IN  IMAGE  SPACE 

CALL  Modulus(ktlpmod,KTLP,n) 

CALL  Modulus(bslpmo(l,BSLP,n) 

c  NORMALIZE  THE  MODULI  FOR  PLOTTING 

CALL  Normalize(ktlpmod,n) 

CALL  Nonnalize(bslpmod,n) 

c  WRITE  THE  MODULUS  TO  A  FILE  FOR  PLOTTING  IF  DESIRED 

20  WRITER, •)’WRnE  MODULUS  TO  A  FILE?  (Y/N)’ 

READ(*,*)rite 
WRITE(*,*)’  ’ 

IF  ((rite.NE.’Y’)j\ND.(rite.NE.’y’).AND. 

+  (rite.NE.’N’).AND,(rite.NE.’n’))  THEN 
WRrTE(*,*)’ERROR,  REENTER.’ 

WRnE(*,*)’  ’ 

GOTO  20 
ENDIF 

IF  ((rite.EQ.’Y’).OR.(rite.EQ.’y’))  THEN 
CALL  Writefile(ktlpmod,bsIpmod,nurapix,lambda,RO,n) 

ENDIF 
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c  DETERMINE  THE  AZIMUTHAL  RMS  PHASE  ERROR 


30  WRITE(*,*)TIND  THE  AZIMUTHAL  RMS  PHASE  ERROR?  (Y/N)’ 
READ(*,*)mean 
WRITE(*,*)’  ’ 

IF  ((mean.NE.’Y’)j\ND.(mean.NE.y)AND. 

+  (mean.NE.’N’)j\ND.(mean.NE.’n’))  THEN 
WRITE(*,*)’ERROR,  REENTER.’ 

WRITE(*,*)’  ’ 

GOTO  30 
ENDIF 

IF  ((mean.EQ.’Y’).OR.(mean.EQ.’y’))  THEN 

c  READ  IN  DATA 

CALL  Readtrufile(TRUDATA,n) 

c  DETERMINE  SQUARE  PHASE  ERROR 

CALL  AMSPE(KTDATA,BSDATA,TRUDATA, 

+  kterror,bseiTor,nyquist,n ; 

c  DETERMINE  AZIMUTHAL  AVERAGE  OF  THE  RMS  PHASE  ERROR 

CALL  AMSPEcalc(kterror,bsenor,rKTeiTor,rBSerror, 

+  nyquist,n) 

c  WRITE  AMSPE  TO  A  FILE  FOR  EACH  CORRELATION  TECHNIQUE 

CALL  WriteiTfile(rKTeiTor,rBSerror,nyquist, lambda, 

+  RO,numpix,iktcount,itccount,n) 

ENDIF 

STOP 

END 
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SUBROUTINE  LIST 


SUBROUTINE  AMSPE(KTDAT^BSDATA,TRUDATA, 

+  kteiTor,bserror,nyqiiist,n) 

c  THIS  S/R  CALCULATES  THE  SQUARE  PHASE  ERROR  OF  THE 

c  RECONSTRUCTED  PHASES 

COMPLEX*  16  KTDATA(n,n),  BSDATA(n,n),  TRUDATA(n,n) 
REAL*8  kterror(n,n),  bseiTor(n,n),  truphase,  ktphase, 

+  bsphase,  pi,  pi2 
pi  =  dacos(-1.0D+00) 
pi2  =  2  •  pi 
n2pl  =  n/2  +  1 
DO  10  i  =  1,  n 
DO  10  j  =  1,  n 
X  =  float(j  -  (n2pl)) 
y  =  float((n2pl)  -  i) 
radius  =  sqrt(x**2.0  +  y**2.0) 

IF  (radius.L£.nyquist)  THEN 
truphase  =  datan2(DIMAG(TRUDATA(ij)), 

+  DREAL(TRUDATA(ij))) 

ktphase  =  datan2(DIMAG(KTDATA(i  j)), 

+  DREAL(KTDATA(y))) 

bsphase  =  datan2(DIMAG(BSDATA(ij)), 

+  DREAL(BSDATA(ij))) 

20  IF  (truphase.GT.pi2)  THEN 

truphase  =  truphase  -  pi2 
GOTO  20 
ENDIF 

30  IF  (truphase.LT.-pi2)  THEN 

truphase  =  truphase  +  pi2 
GOTO  30 
ENDIF 

40  IF  (ktphase.GT.pi2)  THEN 

ktphase  =  ktphase  •  pi2 
GOTO  40 
ENDIF 

50  IF  (ktphase.LT-pi2)  THEN 

ktphase  =  ktphase  +  pi2 
GOTO  50 
ENDIF 
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60  IF  (bsphase.GT.pi2)  THEN 
bsphase  =  bsphase  -  pi2 
GOTO  60 
ENDIF 

70  IF  (bsphase.LT.-pi2)  THEN 

bsphase  =  bsphase  +  pi2 
GOTO  70 
ENDIF 

kteiTor(ij)  =  truphase  -  ktphase 

80  IF  (kterror(ij).GT.pi2)  THEN 

kterror(ij)  =  kterror(ij)  -  pi2 
GOTO  80 
ENDIF 

90  IF  (kterror(ij).LT.-pi2)  THEN 

kterror(ij)  =  kterror(ij)  +  pi2 
GOTO  90 
ENDIF 

IF  (kterror(ij).GT.pi) 

+  kteiTor(ij)  =  IrteiTor(ij)  -  pi2 

IF  (kterror(iJ).LT.-pi) 

+  kteiTor(ij)  =  Irterror(ij)  +  pi2 

kterror(ij)  =  (kteiTor(iJ))**2.0 

bseiTor(ij)  =  truphase  -  bsphase 

100  IF  (bserror(iJ).GT.pi2)  THEN 

bserror(ij)  =  bserror(ij)  -  pi2 
GOTO  100 
ENDIF 

no  IF  (bserror(ij).LT.-pi2)  THEN 

bserror(ij)  =  bserror(ij)  +  pi2 
GOTO  110 
ENDIF 

IF  (bserror(ij).GT.pi) 

+  bserror(ij)  =  bseiTor(ij)  -  pi2 

IF  (bserror(iJ).LT.-pi) 

+  bserror(ij)  =  bseiTor(ij)  +  pi2 

bserror(ij)  =  (bserror(ij))**2.0 
ENDIF 

10  CONTINUE 
RETURN 
END 
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SUBROUTINE  AMSPEcalc(kteiTor,bserror,rKTerror, 

+  rBSerror^yquist,n) 

c  THIS  S/R  CALCULATES  THE  RMS  PHASE  ERROR  AS  A  FUNCTION  OF 

c  RADIUS 

REAL*8  kteiTor(n,n),  bserror(n,n),  rKTerror(ii/2), 

+  rBSerror(n/2) 

INTEGER  r 
n2pl  =  11/2+1 
DO  10  r  =  0,  i^quist 
nerr  =  0 
DO  20  i  =  1,  n 
DO  20  j  =  1,  n 
X  =  float(j  -  (n2pl)) 
y  =  £loat((n2pl)  -  i) 
radius  =  sqrt(x**2.0  +  y**2.0) 

IF  ((radius.GE.float(r))AND. 

+  (radius.LT.float(r+l)))  THEN 
nerr  =  nerr  +  1 

rKTerror(r)  =  rKTerror(r)  +  kterror(ij) 
rBSerror(r)  =  rBSerror(r)  +  bserror(ij) 

ENDIF 

20  CONTINUE 

rKTerror(r)  =  dsqrt(rKTerror(r)/nerr) 
rBSerror(r)  =  dsqrt(rBSerror(r)^err) 

10  CONTINUE 
RETURN 
END 


SUBROUTINE  Modulus(mod,DATA,n) 

c  THIS  S/R  DETERMINES  THE  MODULUS  OF  A  COMPLEX  NUMBER 

COMPLEX*16  DATA(n,n) 

REAL*8  mod(n,n) 

DO  10  i  =  1,  n 
DO  10  j  =  1,  n 
mod(ij)  =  ABS(DATA(iJ)) 

10  CONTINUE 
RETURN 
END 
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SUBROUTINE  Normalize(data,n) 

c  THIS  S/R  NORMALIZES  THE  OUTPUT  DATA  TO  1.0. 

REAL*8  data(n,n),  maxval 
maxval  =  O.OD+00 
DO  10  i  =  1,  n 
DO  10  j  =  1,  n 

IF  (data(iJ).GT.maxval)  maxval  =  data(ij) 

10  CONTINUE 
DO  20  i  =  1,  n 
DO20j  =  1,  n 
data(ij)  =  data(ij)/maxval 
20  CONTINUE 
RETURN 
END 


SUBROUTINE  Readfile(KTDATA,BSDATA,n) 
c  THIS  S/R  READS  THE  RECONSTRUCTED  DATA  FROM  A  FILE 


10 


COMPLEX*16  KTDATA(n,n),  BSDATA(n,n) 

CHARACTER*  16  ktfile,  bsfile 
WRTTEC*,*)’  ’ 

WRrTE(*,*)’ENTER  INPUT  KT  RECON  FILE  NAME  (16  CHAR’ 
WRrTE(*,*)’MAX):’ 

READ  (*,30)  ktffle 
WRrTE(*,*)’  ’ 

WRITER, *)’ENTER  INPUT  TC  RECON  FILE  NAME  (16  CHAR’ 
WRrTE(*,*)’MAX):’ 

READ  (*,30)  bsffle 
WRITE(*,*)’  ’ 

OPEN(UNrT=40,FILE=ktfile,STATUS=’OLD’) 

OPEN(UNrT=50,FILE=bsfile,STATUS=’OLD’) 

DO  10  i  =  1,  n 


DO  10  j  =  1,  n 
READ(40,*)  KTDATA(io) 
CONTINUE 
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DO  20  i  =  1,  n 
DO20j  =  l,n 
READ(50,*)  BSDATA(ij) 
20  CONTINUE 
30  FORMAT(A16) 

RETURN 

END 


SUBROUTINE  Readsiir(ktsnr,bssnr,ktr,bsr4iyquist,n) 

c  THIS  S/R  READS  THE  SNR  DATA  FROM  A  FILE 

REAL*8  ktsnr(ii/2),  bssnr(n/2) 

REAL  ktr(n/2),  bsr(n/2) 

CHARACTER*  16  ktsiirfile,  bssnrfile 

WRITE(*,*)’ENTER  KT  SNR  FILE  NAME  (16  CHAR  MAX):’ 
READ  (*.30)  ktsnrfile 
WRnE(*,*)’  ’ 

WRITE(*,*)’ENTER  TC  SNR  FILE  NAME  (16  CHAR  MAX):’ 
READ  (*,30)  bssnrfile 
WRnE(*,*)’  ’ 

OPEN(UNIT=404TLE=ktsnrfile,STATUS=’OLD’) 

OPEN(UNIT=50,FILE«bssnrfilc,STATUS*’OLD’) 

DO  10  i  =  1,  nyquist 
READ(40,*)  ktr(i),  ktsnr(i) 

10  CONTINUE 

DO  20  i  =  1,  nyquist 
READ(50,*)  bsr(i),  bssnr(i) 

20  CONTINUE 
30  FORMAT(A16) 

RETURN 

END 


SUBROUTINE  Readtrufile(TRUDATAn) 

c  THIS  S/R  READS  THE  TRUTH  DATA  FROM  A  FILE  < 

COMPLEX*  16  TRUDATA(n,n) 

CHARACTER*  16  file 

WRITE(*,*)’ENTER  INPUT  TRUTH  DATA  FILE  NAME  (16  CHAR’ 
WRITE(*,*)’MAX):’ 
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READ  (*,20)  file 
WRITE(*,*)’  ’ 

OPEN(UIOT=30JFILE=file,STATUS=’OLD’) 
DO  10  i  =  1,  n 
DO  10  j  *=  l,n 
READ(30,*)  TRUDATA(ij) 

10  CONTINUE 
20  FORMAT(A16) 

RETURN 

END 


SUBROUTINE  Truncate(DATA,LPDATA,trunc,lambda, 

+  R0,numpix4i) 

c  THIS  S/R  FILTERS  THE  OBJECT  SPECTRUM  ARRAY  BY 
c  TRUNCATING  THE  ARRAY  AT  A  RADIUS  SET  BY  THE  SNR  VALUE 
c  GREATER  THAN  1 

COMPLEX*  16  DATA(n,n),  LPDATA(n,n) 

REAL  lambda 
n2  =  11/2 
n2pl  =  n2  +  1 
pi  =  acos(-1.0E+00) 

DO  10  i  =  1,  n 
DO  10  j  =  l,n 
X  =  float(j-n2pl) 
y  =  float(i-n2pl) 

rad  =  sqrt(x**2  +  y**2)  *  numpix  *  (lambda/RO)  * 

+  (180.0/pi)  *  3600.0 

IF  (rad.LE.trunc)  THEN 
LPDATA(ij)  =  DATA(ij) 

ELSE 

LPDATA(ij)  =  (0.0, 0.0) 

ENDIF 

10  CONTINUE 
RETURN 
END 
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SUBROUTINE  Truncval(snr,r,trunc,nyquist,icount,n) 


c  THIS  S/R  DETERMINES  THE  TRUNCATION  VALUE  BASED  ON  THE 
c  SNR.  THE  VALUE  IS  BASED  ON  THE  POINT  AT  WHICH  THE  SNR  IS 
c  UNITY  OR  GREATER. 

REAL*8  snr(n/2) 

REALr(n/2) 

j  =  1 

DO  10  i  =  1,  nyquist 

IF  ((snr(i).GE.1.0D+00).AND.O.EQ.i))  THEN 
trunc  =  r(i) 
icount  =  j 
j  =  j  +  1 
ENDIF 

10  CONTINUE 
RETURN 
END 


SUBROUTINE  Writefile(ktlpdata,bslpdata,numpix, 
+  lambda,R0,n) 

c  THIS  S/R  WRITES  THE  DATA  TO  A  FILE 


REAL*8  ktlpdata(n,n),  bslpdata(n,n) 

REAL  conv,  lambda,  RO,  x,  y 
CHARACTER*  16  ktlpfile,  bslpflle 
pi  =  acos(-1.0E+00) 

WRITE(*,*)’ENTER  OUTPUT  KT  FILE  NAME  (16  CHAR  MAX):’ 
READ  (*,50)  ktlpfile 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  OUTPUT  TC  FILE  NAME  (16  CHAR  MAX):’ 
READ  (*,50)  bslpfile 
WRITE(*,*)’  ’ 

OPEN(UNIT= 30,FILE= ktlpfile,STATUS= ’NEW’) 
OPEN(UNIT=40,FILE=bslpfile,STATUS=’NEW’) 

DO  10  i  =  1,  n 
DO  10  j  -  1,  n 

conv  =  float(numpix)  *  (lambda/RO)  * 

+  (180.0/pi)  *  3600.0 

X  =  float(j  -  (n/z+l))  *  conv/float(n) 
y  =  float(i  -  (n/2+ 1))  *  conv/float(n) 
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WRITE(30,60)  X,  y,  ktlpdata(y) 

10  CX)NTINIJE 
DO  20  i  =  1,  n 
DO  20  j  =  1,  n 

conv  =  float(numpix)  *  (lambda/RO)  * 
+  (180.0/pi)  •  3600.0 

X  =  float(}  -  (n^+1))  •  conv/float(n) 
y  =  float(i  -  (n/2+1))  •  conv/float(n) 
WRITE(40,60)  X,  y,  bslpdata(ij) 

20  CONTINUE 
50  FORMAT(A16) 

60  FORMAT(F7.4,2x,F7.4,2x,F7.4) 

RETURN 
END 


SUBROUTINE  Writerrfile(rKTerror,rBSeiTor,nyquist, 

+  lambda,RO,numpix,iktcount, 

+  itccount,n) 

c  THIS  S/R  WRITES  THE  ERROR  DATA  TO  A  FILE 

REAL*8  rKTerror(n/2),  rBSeiTor(n^) 

REAL  X,  y,  lambda 
CHARACTER*  16  ktfile,  bsfile 
INTEGER  r 
pi  =  acos(-1.0+00) 

WRITE(*,*)’ENTER  OUTPUT  KT  ERROR  FILE  NAME  (16  CHAR’ 
WRITE(*,*)’MAX);’ 

READ  (*,30)  ktfile 
WRITE(*,*)’  ’ 

WRITE(*,*)’ENTER  OUTPUT  TC  ERROR  FILE  NAME  (16  CHAR’ 
WRITE(*,*)’MAX):’ 

READ  (*,30)  bsfile 
WRITE(*,*)’  ’ 

OPEN(UNIT=  l,FILE=ktfile,STATUS=’NEW’) 
OPEN(UNIT=2,FILE=bsfile,STATUS= ’NEW’) 

DO  10  r  =  1,  iktcount 

X  =  r  *  numpix  *  (lambda/RO)  *  (180.0/pi)  *  3600.0 
WRrrE(l,*)  X,  rKTerror(r) 

10  CONTINUE 
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DO  20  r  =  1,  itccount 

X  =  r  *  numpix  *  (lambda/RO)  *  (180.0/pi)  *  3600.0 
WRITE(2,*)  X,  rBSerror(r) 

20  CONTINUE 
30  FORMAT(A16) 

RETURN 

END 
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